Abundance of Deer in Victoria

Determining abundance of deer species in Victoria using camera trap data

Methods used to estimate deer density

Camera detection

The density of deer at a given site can be estimated using camera-trap distance sampling (CTDS). CTDS is a modified form of distance-sampling that allows us to infer the probability a given individual will be detected within the survey area (area in front of the camera). This detection probability is a function of the distance of the individual from the camera, whereby individuals entering the camera field of view further from the camera have a lower detection probability up to a given truncation distance from the camera where detection probability is near-zero.

An underlying assumption about CTDS is that the probability a deer will be available for detection at any given point location within the camera field of view is equal. Under this assumption, for a point transect, we take into account the total area for each distance-bin area, which increases at further distance bins. However, in this study we implement a novel method that considers group size of the detected species in the availability calculations. For larger groups, CTDS should account for the availability of the closest individual rather than the availability of all individuals. When group size increases it is more likely the the closest individual from the group is closer to the camera. If we assume that the triggering of the camera is dependent upon the closest individual than we must adjust our estimated availability to account for variable group sizes. If we do not adjust for group size and only use the distance to the closest individual for our DS models then we will likely underestimate the detection rate as distances will be smaller than reality. Alternatively, if we record distances to multiple individuals in the same photo and take an average or model them independently we would likely overestimate detection probability because individuals at further distances are recorded because a closer individual has triggered the camera trap. For more information on how group size is likely to effect CTDS abundance estimates a simulation study has been written here: https://justincally.github.io/blog/posts/2022-11-10-ctds-for-groups/

In this study we investigate two possible detection functions that may explain how detection rates ‘fall off’ over the distances from the camera. Firstly a half-normal detection function given by:

\[p = exp(-y^2/2\sigma^2) \] Hereby \(p\) is the probability of detection, \(y\) is the distance from the camera (midpoint of the detection bin), and \(\sigma\) is the shape parameter. Alternatively, we compare the fit of this detection function to a hazard-rate function, where detection probability is given by:

\[p = 1 - exp(-(y/\sigma)^{-\theta})\] In these hazard-rate models, an additional parameter (\(\theta\)) is estimated, which is scale parameter. Given that detection rate may not necessarily be uniform across sites we model the magnitude of the shape parameter across sites (i) as:

\[log(\sigma_i) = \alpha + \beta_{det} \times HU_i\] Where \(HU_i\) is the amount of herbaceous understorey measured surrounding the camera. \(\beta_{det}\) is the coefficient for the effect of herbaceous understorey on detection and \(\alpha\) is the intercept.

Using these estimates of detection probability in different bins we can then estimate the average detection probability of an individual at a given camera station. This can then be used to account for imperfect detection when fitting observed counts from the camera trap to parameters estimating the true abundance at a site.

Modelling density

Modelling abundnace using camera trap data

The density at a site is modelled as a function of the observed counts, relative frequency of group sizes, distance-sampling detection probability, survey effort (area in front of camera x time/snapshot moments the camera is deployed for), proportion of time within a 24-hour cycle that deer were active (\(activity\)). This parameter was estimated using the activity R package and included in the model using an informative beta-distribution prior. For a given species the number of detection’s of groups (1 or more deer) during the deployment is estimated with:

\[C_{n,j} = NegBinomial2(exp(log(\lambda_{n}) + log(p_{n,j}) + log(activity) + log(\epsilon gs_{n,j}) + log(\bar p)), \phi_n) \times survey\ area_n\]
Where \(C_{n,j}\) is the observed count of group size (\(j\)) at a site (\(n\)). \(\lambda_{n}\) is the true density at a site. In this model \(X_n\) are the values of the fixed-effect covariates, \(\beta\) are the fixed effect coefficients and \(\epsilon bioregion_n\) is the bioregion random effect drawn from a standard normal distribution:

\[\lambda_{n} = X_n \times \beta + \epsilon bioregion_n\]

The \(\epsilon gs_{n,j}\) is the proportional variation in group size ar a site, where \(\sum \epsilon gs_{n,j \in [1,max(group\ size)]} = 1\). Assuming that the rate of group sizes can vary between sites, we model the variation in group size between sites with a group-size level intercept (\(\zeta_j\)) and site-group-size level random effect (\(\epsilon site_{n,j}\)):

\[\epsilon psi_{nj} = exp(\zeta_j + \epsilon site_{n,j}) \] The proportional group size at a given site and for a given group size is thus given by:

\[\epsilon gs_{nj} = \frac{\epsilon psi_{nj}}{sum(\epsilon psi_{n})}\]

We opt for a negative binomial model to account for over-dispersion in counts/abundance. The parameterization of the negative binomial model implemented here (NegBinomial2) contains a mean as expressed by the log-sum of the first five terms. It also contains a dispersion parameter (\(\phi_n\)), which is the over-dispersion at the site. \(\phi_n\) is modelled as a function of bioregion:

\[\phi_n = exp(\chi mu_{bioregion_n} + \chi sd_{bioregion_n} \times \chi raw_{bioregion_n})\]

Supplementing the camera trap data with transect-level data

Our surveys involved the deployment of cameras at a site for usually 6-8 weeks as well as 150m transect searches that were walked up and back by a single individual for two transects, and one transect which was walked once by two observers. This gives us an approximate survey area at each site of 150m x 2 x 3 = 900m. Deer sign was recorded on these transects with pellets, footprints, rubbings, and wallows recorded as either present or absent.

This data provides supplemental presence-absence data that can be integrated into our model to help inform likely densities at sites where cameras did not record observations but one or more deer signs were detected along the transects.

The foundation for this integration is a Royle-Nichols model (Royle and Nichols 2003) that models the frequency of detection as a mixture of detection probability and abundance. Given that we also collect absolute abundance estimates via the camera trap data, relative abundance measures from this output can be transformed to absolute abundance and thus provide an avenue to integrate camera and transect-based data into a single model estimating abundance. The model can be expressed as:

\[N_n \sim NegBinomial2(\lambda_n, \bar\phi)\] Where \(\bar\phi\) is the over-dispersion term of the Royle-Nichols component of the model. Conditional detection probability is expressed based on the probability of detection from a given visit/survey type (\(j\)) (e.g. camera, pellet transect, footprint transect), and the number of individuals available to be detected at the site (\(N_n\)):

\[p_{nj} = 1 - (1 - r_{nj}) ^ {N_n}\] The model for this detection probability is then:

\[y_{nj} \sim Bernoulli(p_{nj})\] This model thus allows for \(\lambda_n\) to be modelled via the camera counts and the results of detections along the transect.

Finally, the Royle-Nichols model also allows us to estimate the rates at which with camera trap did not detect individuals at a site when they were present. Specifically we can use the average probability of detection of at least one individual at the camera trap, given individuals have expected abundance lambda (\(\bar p\)).

Setup

Load in relevant packages for analysis, additionally, connect to the database. Camera trap and site data is stored on the database.

Show code

library(cmdstanr)
library(dplyr)
library(sf)
library(bayesplot)
library(loo)
library(gt)
library(terra)
library(caret)
library(tidyterra)
library(tidyr)
library(VicmapR)
library(kableExtra)
library(brms)
library(Distance)
options(brms.backend = "cmdstanr")
check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)
register_knitr_engine(override = FALSE)
options(mc.cores=8)

#### Database connection ####
con <- weda::weda_connect(password = keyring::key_get(service = "ari-dev-weda-psql-01",
                                                username = "psql_user"), username = "psql_user")

Custom Functions

Additional functions used in the data preparation, modelling and analysis are available in the /functions directory.

Show code

# Source internal functions
sapply(list.files("functions", full.names = T), source, verbose = F)

Prepare Data

Wrangle and format data for the STAN models for the various species.

Scope of models

Outline which species should be modeled, and which projects to source data from.

Show code

# Species to run model for.
deer_species_all <- c("Cervus unicolor", "Dama dama", "Cervus elaphus", "Axis porcinus")
species_names <- c("Sambar Deer", "Fallow Deer", "Red Deer", "Hog Deer")
# Projects to select.
project_short_name <- c("hog_deer_2023", "StatewideDeer")
# Buffer for data extraction.
spatial_buffer <- 1000
# Covariate Rasters
raster_files <- "/Volumes/DeerVic\ Photos/Processed_Rasters"
# raster_files <- "data/prediction_raster"
prediction_raster <- "data/prediction_raster/statewide_raster.tif"
# For the integrated model we place limits on the maximum density of deer to integrate over. In cases where there are no detections on the camera this is limited to 15 per km2. In cases where deer were detected on the camera this is expanded to 50. We believe this is sufficiently high. Very high values of these will be less efficient. 
n_max_no_det <- 15
n_max_det <- 30

Camera locations

Download the camera locations from the database, this table outlines the locations and the deployment history of the cameras.

Show code

cams_curated <- tbl(con, dbplyr::in_schema("camtrap", "curated_camtrap_operation")) %>%
  dplyr::filter(ProjectShortName %in% !!project_short_name) %>% # Only retrieve cam records for hog deer 2023
  dplyr::collect() %>%
  dplyr::arrange(SiteID) %>%
  sf::st_as_sf(., coords = c("Longitude", "Latitude"), crs = 4283) 

n_site <- nrow(cams_curated)

Formulas for detection and abundance

Here we outline formulas to be used in the models. The formulas account for the various fixed-effect parameters.

Show code

#### Model formulas ####

#### Transect Formula: Survey Only ####
transect_formula <- ~Survey

#### Abundance Formula Options #### 
ab_formula_1 <- ~ scale(sqrt(BareSoil)) + scale(sqrt(Nitrogen)) + scale(log(PastureDistance)) + scale(BIO15) + scale(sqrt(ForestEdge)) + 
  scale(sqrt(MRVBF)) + scale(sqrt(SLOPE)) + scale(sqrt(NonPhotosyntheticVeg))

ab_formula_2 <- ~ scale(sqrt(BareSoil)) + scale(sqrt(Nitrogen)) + scale(log(PastureDistance)) + scale(BIO15) + scale(sqrt(ForestEdge)) + scale(sqrt(MRVBF)) 

ab_formula_3 <- ~ scale(sqrt(BareSoil)) + scale(sqrt(Nitrogen)) + scale(log(PastureDistance)) + scale(BIO15) + scale(sqrt(ForestEdge)) 

ab_formula_4 <- ~ 1

#### Detection Formula: Distance-sampling ####
det_formula <- ~ scale(HerbaceousUnderstoryCover) # average detection across all sites 

Create model data

Using the prepare_model_data() function we generate the data for the STAN model. This function will:
1. Download data from the database 2. Format that data to match the distance-sampling bins
3. Organised the counts into sites and group sizes
4. Generate model matrix of the various sub-models (distance-sampling, abundance and transect) 5. Generate data for the prediction process 6. Generate data for the random effect (bioregion)
7. Generate data for regional predictions (indexing based on DEECA regions)

Show code

if(!file.exists("data/multispecies_data.rds")) {
  
crown_land <- readRDS("data/crown_land.rds")
  
multispecies_data <- prepare_model_data_multispecies(species = deer_species_all,
                                 projects = project_short_name,
                     buffer = spatial_buffer,
                     detection_formula = det_formula,
                     abundance_formula = ab_formula_1,
                     transect_formula = transect_formula,
                     con = con,
                     raster_dir = raster_files,
                     prediction_raster = prediction_raster,
                     n_max_no_det = n_max_no_det,
                     n_max_det = n_max_det,
                     crown_land = crown_land,
                     evaltransects = TRUE, 
                     filter_behaviour = TRUE)

multispecies_data$keyfun <- 1 # 0 = HN, 1 = HZ
multispecies_data$raw_data[is.na(multispecies_data$raw_data)] <- 0
multispecies_data$transects[is.na(multispecies_data$transects)] <- 0

evc_groups <- readRDS("data/evc_groupnames.rds")

multispecies_data <- c(multispecies_data, evc_groups)

saveRDS(multispecies_data, "data/multispecies_data.rds")
} else {
  multispecies_data <- readRDS("data/multispecies_data.rds")
}

Model Execution

MCMC settings

Below we list the MCMC setting for our model. We run models on eight parallel chains for 400 iterations each (200 warm-up and 200 sampling). These setting provide us with 1,600 posterior draws (8 x 200).

Show code

# STAN settings
ni <- 200 # sampling iterations
nw <- 200 # warm-up iterations 
nc <- 8 # number of chains

Read in Models

These are the models used in the analysis. The first is an integrated model that requires transect and camera data. The second is a count-only model that just requires the camera data.

Show code

functions {

    /* Half-normal function
   * Args:
   *   sigma: sigma
   *   midpoints: midpoints
   * Returns:
   *   detection probability
   */
    vector halfnorm(real sigma, vector midpoints) {
      int bins = rows(midpoints);
      vector[bins] p_raw; // detection probability

      p_raw = exp(-(midpoints)^2/(2*sigma^2));
      return p_raw;
    }

      /* Hazard function
   * Args:
   *   sigma: sigma
   *   theta: theta
   *   midpoints: midpoints
   * Returns:
   *   detection probability
   */
    vector hazard(real sigma, real theta, vector midpoints) {
      int bins = rows(midpoints);
      vector[bins] p_raw; // detection probability

      p_raw = 1 - exp(-(midpoints/sigma)^(-theta));
      return p_raw;
    }

  vector prob_dist(real sigma, real theta, int keyfun, vector midpoints){
  int bins = rows(midpoints);
  vector[bins] out; // detection probability

  if(keyfun == 0){
    out = halfnorm(sigma, midpoints);
  } else if(keyfun == 1){
    out = hazard(sigma, theta, midpoints);
  }
  return out;
  }
}
data {
  int<lower=0> N;                      // number of observations
  int<lower=0> S;                      // number of species
  real delta;                          // bin width
  int<lower=1> n_site;                 // site
  int<lower=1> n_distance_bins;        // distance bins
  int<lower=1> n_gs;                   // number of group sizes
  vector[n_gs] gs;                     // group sizes
  vector[n_distance_bins] midpts;      // midpt of each interval
  real<lower=1> max_distance;         // truncation distance
  int<lower=1> max_int_dist;          // max distance as an integer
  real<lower=0> theta_frac;           // fraction of camera view
  array[n_site] int effort;           // effort
  array[n_site, n_gs, S] int n_obs;      //number of observations
  array[n_site, n_distance_bins, n_gs] int y; // observations matrix

  // summary of whether species is known to be present at each site
  int<lower = 0, upper = 1> any_seen[S, n_site];

  // number of surveys at each site
  int<lower = 0> n_survey[n_site];

  // availability prior information
  real<lower=0> bshape;               // availability shape
  real<lower=0> bscale;               // availability scale

  // camera trap detection parameters
  int<lower=0> det_ncb;                 // number of covariates for detection model
  matrix[n_site, det_ncb] det_model_matrix; // detection model matrix
  array[n_gs, n_distance_bins] real pa ; // averaged availability for multiple individuals

  // Abundance/occupancy model matrix
  int<lower = 1> m_psi;
  matrix[n_site, m_psi] X_psi;
  // negbinom scale
  // real reciprocal_phi_scale;

  //transect level information
  int<lower=1> trans;                  // total number of transects across all sites for all methods
  array[S, trans] int<lower = 0, upper = 1> y2; // transect based binomial detections
  int<lower = 0, upper = trans> start_idx[n_site];
  int<lower = 0, upper = trans> end_idx[n_site];
  int<lower=1> trans_det_ncb;           // number of covariates for transect detection model
  matrix[trans, trans_det_ncb] trans_pred_matrix; // transect detection model matrix
  int<lower=1>  n_max[n_site, S]; // max for poisson RN

  // Prediction data
  int<lower=1> npc;                 // number of prediction grid cells
  matrix[npc, m_psi] X_pred_psi; // pred matrix
  vector[npc] prop_pred; //offset
  // bioregion RE
  int<lower=1> np_bioreg;
  int<lower=1> site_bioreg[n_site];
  int<lower=1> pred_bioreg[npc];
    // region data
  int<lower=1> np_reg;
  int<lower=1> site_reg[n_site];
  int<lower=1> pred_reg[npc];

  // key function, 0 = halfnorm
  int keyfun;
}

transformed data {
  // vector[n_distance_bins] pi; // availability for point
  vector[n_site] survey_area;
  array[S, n_site] real cam_seen;
    for(i in 1:n_site) {
      survey_area[i] = theta_frac * effort[i] * pi() * (max_distance/1000)^2;
      for(s in 1:S) {
      cam_seen[s,i] = sum(n_obs[i,,s]);
      }
    }

}

parameters {
 // abundance parameters
  //array[n_site] simplex[n_gs] eps_ngs; // random variation in group size
  array[S] vector[m_psi] beta_psi;
  vector[det_ncb] beta_det;
  real log_theta;
  // transect detection parameters
  vector[trans_det_ncb] beta_trans_det;
  // temporal availability parameters
  real<lower=0, upper=1> activ;
  // bioregion RE
  array[S] real<lower=0> bioregion_sd;
  array[S] vector[np_bioreg] bioregion_raw;
  // eps group size params
  array[S] vector[n_gs] zeta;
  array[S] matrix[n_site, n_gs] eps_raw;
  array[S] real<lower=0> grp_sd;
  // od
  array[S] real od_mu;
  array[S] real<lower=0> od_sd;
  array[S] vector[np_bioreg] od_raw;
  real odRNmu;
}

transformed parameters {
  // random effects
  array[S] vector[np_bioreg] eps_bioregion; // bioregion random effect
  // distance parameters
  array[n_site] real log_sigma;
  array[n_site] real sigma;
  array[n_site] vector[n_distance_bins] p_raw; // detection probability
  array[n_site, n_distance_bins, n_gs] real p_raw_scale; // detection probability for point independence model and accounting for availability
  array[n_site, n_distance_bins, n_gs] real<upper=0> log_p_raw; // log detection probability accounting for availability
  array[n_site, n_gs] real log_p; // log probability with summed exponential for the multinomial logit model
  array[n_site, n_gs] real<lower=0,upper=1> p; // log probability with summed exponential for the multinomial logit model
  // abundance parameters
  array[S, n_site, n_gs] real<lower=0> lambda;
  // activity parameters
  real log_activ = log(activ);
  vector[trans] logit_trans_p = trans_pred_matrix * beta_trans_det; // observation process model
  // lp_site for RN model
  array[S, n_site] real lp_site;
  vector<lower=0,upper=1>[trans] r = inv_logit(logit_trans_p);
  array[S] vector[n_site] log_lambda_psi;
  // eps group size
  array[S, n_site] simplex[n_gs] eps_ngs; // random variation in group size
  array[S, n_site] vector[n_gs] epsi;
  array[S] matrix[n_site, n_gs] eps_site;
  real<lower=0> theta = exp(log_theta);
  array[S] vector[np_bioreg] od; // bioregion random effect
  real odRN = exp(odRNmu);  // bioregion random effect
  array[S, n_site] real<lower=0, upper=1> pbar;

  for(s in 1:S) {
    for(b in 1:np_bioreg) {
    eps_bioregion[s,b] = bioregion_sd[s] * bioregion_raw[s,b];
    od[s,b] = exp(od_mu[s] + od_sd[s] * od_raw[s,b]);
  }
  }



for(n in 1:n_site) {
  log_sigma[n] = det_model_matrix[n,] * beta_det;
  sigma[n] = exp(log_sigma[n]);
  p_raw[n] = prob_dist(sigma[n], theta, keyfun, midpts);
  for (i in 1:n_distance_bins) {
      // assuming a half-normal detection fn from line
      // p_raw[n,i] = exp(-(midpts[i])^2/(2*sigma[n]^2)); // half-normal function (pg 419 of AHM)
      // hazard function
      // p_raw[n,i] = 1 - exp(-(midpts[i]/theta)^(-sigma[n])); // hazard function (pg 507 of AHM)
    for(j in 1:n_gs) {
      p_raw_scale[n,i,j] = p_raw[n,i]*pa[j,i]; //  pr(animal occurs and is detected in bin i)
      log_p_raw[n,i,j] = log(p_raw_scale[n,i,j]);
      }
  }

  for(s in 1:S) {
  vector[n_max[n,s]+1] Nlik;
  vector[n_max[n,s]+1] gN;
  vector[n_max[n,s]+1] pcond;
  log_lambda_psi[s,n] = X_psi[n,] * beta_psi[s] + eps_bioregion[s,site_bioreg[n]];

  for(j in 1:n_gs) {
    eps_site[s, n,j] = grp_sd[s] * eps_raw[s,n,j];
    epsi[s,n,j] = exp(zeta[s,j] + eps_site[s,n,j]);
  }
  eps_ngs[s,n,] = epsi[s,n,]/sum(epsi[s,n,]);

  // p-bar

  for(k in 1:(n_max[n,s]+1))
    Nlik[k] = exp(neg_binomial_2_log_lpmf(k-1 | log_lambda_psi[s,n], odRN));

  gN = Nlik/sum(Nlik);

  for(k in 1:(n_max[n,s]+1))
    pcond[k] = (1 - (1 - r[start_idx[n]])^(k-1)) * gN[k];

  pbar[s,n] = sum(pcond);

  }

  for(j in 1:n_gs) {
    log_p[n,j] = log_sum_exp(log_p_raw[n,,j]);
    p[n,j] = exp(log_p[n,j]);
    // model site abundance
    for(s in 1:S) {
    lambda[s,n,j] = exp(log_lambda_psi[s,n] + log_p[n,j] + log_activ + log(eps_ngs[s,n,j]) + log(pbar[s,n])) .* survey_area[n];
    }
  }

// Royle-Nichols implementation in STAN (looping over possible discrete values of N)
// https://discourse.mc-stan.org/t/royle-and-nichols/14150
// https://discourse.mc-stan.org/t/identifiability-across-levels-in-occupancy-model/5340/2
    for(s in 1:S) {
if (n_survey[n] > 0) {
  vector[n_max[n,s]] lp;
    if(any_seen[s,n] == 0){ // not seen
      lp[1] = log_sum_exp(neg_binomial_2_log_lpmf(0 | log_lambda_psi[s,n], odRN), neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]));
    } else {
      lp[1] = neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]);
    }
    for (k in 2:n_max[n,s]){
      lp[k] = neg_binomial_2_log_lpmf(k | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | 1-(1-r[start_idx[n]:end_idx[n]])^k);
    }
    lp_site[s,n] = log_sum_exp(lp);
    } else {
    lp_site[s, n] = 0;
    }
  }
    }
}

model {

  for(s in 1:S) {
    beta_psi[s] ~ normal(0, 3); // prior for intercept in poisson model
    bioregion_sd[s] ~ normal(0,2);
    bioregion_raw[s,] ~ normal(0,1);
    to_vector(eps_raw[s,,]) ~ std_normal();
    grp_sd[s] ~ normal(0, 1);
    zeta[s,] ~ normal(0, 2);
    // od
    od_sd[s] ~ normal(0,1);
    od_mu[s] ~ normal(0,1);
    od_raw[s,] ~ normal(0,1);
  }
  odRNmu ~ normal(0,1);
  beta_trans_det ~ normal(0, 2); // beta for transect detection
  beta_det ~ normal(0, 4); // prior for sigma
  activ ~ beta(bshape, bscale);  //informative prior
  log_theta ~ normal(0,2); // prior for theta

  for(n in 1:n_site) {
  for(j in 1:n_gs) {
  for(s in 1:S) {
  target += neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s, site_bioreg[n]]);
        }
  y[n,,j] ~ multinomial_logit(to_vector(log_p_raw[n,,j]));
  }
  for(s in 1:S) {
  target += lp_site[s, n];
  }
}
}

generated quantities {
  array[S, n_site,n_gs] real n_obs_pred;
  array[S, n_site, n_gs] real n_obs_true;
  array[S, n_site] real N_site;
  array[S, n_site] real N_site_pred;
  array[n_site, max_int_dist+1] real DetCurve;
  array[n_site, n_gs] real log_lik1;
  array[S, n_site, n_gs] real log_lik2;
  array[n_site, n_gs] real log_lik2_site;
  array[S, n_site] real log_lik2_species;
  array[n_site] real log_lik;
  array[n_site] real log_lik_det;
  //array[S, n_site] real Site_lambda;
  real av_gs[S];
  array[S] simplex[n_gs] eps_gs_ave;
  array[S, npc] real pred;
  //array[S, npc] real pred_trunc;
  array[S, np_reg] real Nhat_reg;
  array[S, np_bioreg] real Nhat_bioreg;
  // array[np_reg] real Nhat_reg_design;
  real Nhat[S];
  //real Nhat_trunc[S];
  real Nhat_sum;
  int trunc_counter[S];
  trunc_counter[S] = 0;

for(s in 1:S) {
  eps_gs_ave[s] = exp(zeta[s])/sum(exp(zeta[s]));
  av_gs[s] = sum(gs .* eps_gs_ave[s]); //average group size
}


for(n in 1:n_site) {
  for(j in 1:n_gs) {
  log_lik1[n,j] = multinomial_logit_lpmf(y[n,,j] | to_vector(log_p_raw[n,,j])); //for loo
  for(s in 1:S) {
  log_lik2[s,n,j] = neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s, site_bioreg[n]]); //for loo
  n_obs_true[s, n, j] = gs[j] * (neg_binomial_2_log_rng(log_lambda_psi[s,n] + log(eps_ngs[s,n,j]), od[s, site_bioreg[n]]));
  n_obs_pred[s, n,j] = gs[j] * neg_binomial_2_rng(lambda[s,n,j], od[s, site_bioreg[n]]);
  }
  log_lik2_site[n, j] = log_sum_exp(log_lik2[,n,j]);
    }
    // get loglik on a site level
    log_lik_det[n] = log_sum_exp(log_lik1[n,]);
    log_lik[n] = log_sum_exp(log_sum_exp(log_lik_det[n],
    log_sum_exp(log_lik2_site[n,])), log_sum_exp(lp_site[,n]));
      for(s in 1:S) {
    //Site_lambda[s,n] = exp(log_lambda_psi[s,n]);
    log_lik2_species[s, n] = log_sum_exp(log_sum_exp(log_lik2[s,n,]), lp_site[s,n]);
    N_site[s,n] = sum(n_obs_true[s,n,]);
    N_site_pred[s,n] = sum(n_obs_pred[s,n,]);
      }

  // loop over distance bins
  if(keyfun == 0) {
  for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
    DetCurve[n, j+1] = exp(-(j+0.5)^2/(2*sigma[n]^2)); // half normal
    }
  } else if(keyfun == 1) {
    for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
    DetCurve[n, j+1] =  1 - exp(-((j+0.5)/sigma[n])^(-theta)); //hazard rate
    }
  }
}

for(s in 1:S) {
for(i in 1:np_reg) Nhat_reg[s,i] = 0;
for(i in 1:np_bioreg) Nhat_bioreg[s,i] = 0;

for(i in 1:npc) {
  pred[s,i] = neg_binomial_2_log_rng(X_pred_psi[i,] * beta_psi[s] + eps_bioregion[s, pred_bioreg[i]], od[s, pred_bioreg[i]]) * prop_pred[i] * av_gs[s]; //offset
  if(pred[s,i] > max(N_site[s,])) {
    //pred_trunc[s,i] = max(N_site[s,]);
    trunc_counter[s] += 1;
  }
  Nhat_reg[s,pred_reg[i]] += pred[s,i];
  Nhat_bioreg[s,pred_bioreg[i]] += pred[s,i];
}
Nhat[s] = sum(pred[s,]);
//Nhat_trunc[s] = sum(pred_trunc[s,]);
}
Nhat_sum = sum(Nhat);
}

Show code

functions {

    /* Half-normal function
   * Args:
   *   sigma: sigma
   *   midpoints: midpoints
   * Returns:
   *   detection probability
   */
    vector halfnorm(real sigma, vector midpoints) {
      int bins = rows(midpoints);
      vector[bins] p_raw; // detection probability

      p_raw = exp(-(midpoints)^2/(2*sigma^2));
      return p_raw;
    }

      /* Hazard function
   * Args:
   *   sigma: sigma
   *   theta: theta
   *   midpoints: midpoints
   * Returns:
   *   detection probability
   */
    vector hazard(real sigma, real theta, vector midpoints) {
      int bins = rows(midpoints);
      vector[bins] p_raw; // detection probability

      p_raw = 1 - exp(-(midpoints/sigma)^(-theta));
      return p_raw;
    }

  vector prob_dist(real sigma, real theta, int keyfun, vector midpoints){
  int bins = rows(midpoints);
  vector[bins] out; // detection probability

  if(keyfun == 0){
    out = halfnorm(sigma, midpoints);
  } else if(keyfun == 1){
    out = hazard(sigma, theta, midpoints);
  }
  return out;
  }
}
data {
  int<lower=0> N;                      // number of observations
  int<lower=0> S;                      // number of species
  real delta;                          // bin width
  int<lower=1> n_site;                 // site
  int<lower=1> n_distance_bins;        // distance bins
  int<lower=1> n_gs;                   // number of group sizes
  vector[n_gs] gs;                     // group sizes
  vector[n_distance_bins] midpts;      // midpt of each interval
  real<lower=1> max_distance;         // truncation distance
  int<lower=1> max_int_dist;          // max distance as an integer
  real<lower=0> theta_frac;           // fraction of camera view
  array[n_site] int effort;           // effort
  array[n_site, n_gs, S] int n_obs;      //number of observations
  array[n_site, n_distance_bins, n_gs] int y; // observations matrix

  // summary of whether species is known to be present at each site
  int<lower = 0, upper = 1> any_seen[S, n_site];

  // number of surveys at each site
  int<lower = 0> n_survey[n_site];

  // availability prior information
  real<lower=0> bshape;               // availability shape
  real<lower=0> bscale;               // availability scale

  // camera trap detection parameters
  int<lower=0> det_ncb;                 // number of covariates for detection model
  matrix[n_site, det_ncb] det_model_matrix; // detection model matrix
  array[n_gs, n_distance_bins] real pa ; // averaged availability for multiple individuals

  // Abundance/occupancy model matrix
  int<lower = 1> m_psi;
  matrix[n_site, m_psi] X_psi;
  // negbinom scale
  // real reciprocal_phi_scale;

  //transect level information
  int<lower=1> trans;                  // total number of transects across all sites for all methods
  array[S, trans] int<lower = 0, upper = 1> y2; // transect based binomial detections
  int<lower = 0, upper = trans> start_idx[n_site];
  int<lower = 0, upper = trans> end_idx[n_site];
  int<lower=1> trans_det_ncb;           // number of covariates for transect detection model
  matrix[trans, trans_det_ncb] trans_pred_matrix; // transect detection model matrix
  int<lower=1>  n_max[n_site, S]; // max for poisson RN

  // Prediction data
  int<lower=1> npc;                 // number of prediction grid cells
  matrix[npc, m_psi] X_pred_psi; // pred matrix
  vector[npc] prop_pred; //offset
  // bioregion RE
  int<lower=1> np_bioreg;
  int<lower=1> site_bioreg[n_site];
  int<lower=1> pred_bioreg[npc];
    // region data
  int<lower=1> np_reg;
  int<lower=1> site_reg[n_site];
  int<lower=1> pred_reg[npc];
  // evc data
  int<lower=1> np_evc;
  int<lower=1> site_evc[n_site];
  int<lower=1> pred_evc[npc];


  // key function, 0 = halfnorm
  int keyfun;
}

transformed data {
  // vector[n_distance_bins] pi; // availability for point
  vector[n_site] survey_area;
  array[S, n_site] real cam_seen;
    for(i in 1:n_site) {
      survey_area[i] = theta_frac * effort[i] * pi() * (max_distance/1000)^2;
      for(s in 1:S) {
      cam_seen[s,i] = sum(n_obs[i,,s]);
      }
    }

}

parameters {
 // abundance parameters
  //array[n_site] simplex[n_gs] eps_ngs; // random variation in group size
  array[S] vector[m_psi] beta_psi;
  vector[det_ncb] beta_det;
  real log_theta;
  // transect detection parameters
  vector[trans_det_ncb] beta_trans_det;
  // temporal availability parameters
  real<lower=0, upper=1> activ;
  // bioregion RE
  array[S] real<lower=0> bioregion_sd;
  array[S] vector[np_bioreg] bioregion_raw;
  // evc RE
  real<lower=0> evc_sd;
  array[S] vector[np_evc] evc_raw;
  // eps group size params
  array[S] vector[n_gs] zeta;
  array[S] matrix[n_site, n_gs] eps_raw;
  array[S] real<lower=0> grp_sd;
  // od
  array[S] real od_mu;
  array[S] real<lower=0> od_sd;
  array[S] vector[np_bioreg] od_raw;
  real odRNmu;
}

transformed parameters {
  // random effects
  array[S] vector[np_bioreg] eps_bioregion; // bioregion random effect
  array[S] vector[np_evc] eps_evc; // bioregion random effect
  // distance parameters
  array[n_site] real log_sigma;
  array[n_site] real sigma;
  array[n_site] vector[n_distance_bins] p_raw; // detection probability
  array[n_site, n_distance_bins, n_gs] real p_raw_scale; // detection probability for point independence model and accounting for availability
  array[n_site, n_distance_bins, n_gs] real<upper=0> log_p_raw; // log detection probability accounting for availability
  array[n_site, n_gs] real log_p; // log probability with summed exponential for the multinomial logit model
  array[n_site, n_gs] real<lower=0,upper=1> p; // log probability with summed exponential for the multinomial logit model
  // abundance parameters
  array[S, n_site, n_gs] real<lower=0> lambda;
  // activity parameters
  real log_activ = log(activ);
  vector[trans] logit_trans_p = trans_pred_matrix * beta_trans_det; // observation process model
  // lp_site for RN model
  array[S, n_site] real lp_site;
  vector<lower=0,upper=1>[trans] r = inv_logit(logit_trans_p);
  array[S] vector[n_site] log_lambda_psi;
  // eps group size
  array[S, n_site] simplex[n_gs] eps_ngs; // random variation in group size
  array[S, n_site] vector[n_gs] epsi;
  array[S] matrix[n_site, n_gs] eps_site;
  real<lower=0> theta = exp(log_theta);
  array[S] vector[np_bioreg] od; // bioregion random effect
  real odRN = exp(odRNmu);  // bioregion random effect
  array[S, n_site] real<lower=0, upper=1> pbar;

  for(s in 1:S) {
    for(b in 1:np_bioreg) {
    eps_bioregion[s,b] = bioregion_sd[s] * bioregion_raw[s,b];
    od[s,b] = exp(od_mu[s] + od_sd[s] * od_raw[s,b]);
  }
  for(k in 1:np_evc) {
  eps_evc[s,k] = evc_sd * evc_raw[s,k];
}
  }

for(n in 1:n_site) {
  log_sigma[n] = det_model_matrix[n,] * beta_det;
  sigma[n] = exp(log_sigma[n]);
  p_raw[n] = prob_dist(sigma[n], theta, keyfun, midpts);
  for (i in 1:n_distance_bins) {
      // assuming a half-normal detection fn from line
      // p_raw[n,i] = exp(-(midpts[i])^2/(2*sigma[n]^2)); // half-normal function (pg 419 of AHM)
      // hazard function
      // p_raw[n,i] = 1 - exp(-(midpts[i]/theta)^(-sigma[n])); // hazard function (pg 507 of AHM)
    for(j in 1:n_gs) {
      p_raw_scale[n,i,j] = p_raw[n,i]*pa[j,i]; //  pr(animal occurs and is detected in bin i)
      log_p_raw[n,i,j] = log(p_raw_scale[n,i,j]);
      }
  }

  for(s in 1:S) {
  vector[n_max[n,s]+1] Nlik;
  vector[n_max[n,s]+1] gN;
  vector[n_max[n,s]+1] pcond;
  log_lambda_psi[s,n] = X_psi[n,] * beta_psi[s] + eps_bioregion[s,site_bioreg[n]] + eps_evc[s,site_evc[n]];

  for(j in 1:n_gs) {
    eps_site[s, n,j] = grp_sd[s] * eps_raw[s,n,j];
    epsi[s,n,j] = exp(zeta[s,j] + eps_site[s,n,j]);
  }

  eps_ngs[s,n,] = epsi[s,n,]/sum(epsi[s,n,]);

  // p-bar
  for(k in 1:(n_max[n,s]+1))
    Nlik[k] = exp(neg_binomial_2_log_lpmf(k-1 | log_lambda_psi[s,n], odRN));

  gN = Nlik/sum(Nlik);

  for(k in 1:(n_max[n,s]+1))
    pcond[k] = (1 - (1 - r[start_idx[n]])^(k-1)) * gN[k];

  pbar[s,n] = sum(pcond);

  }

  for(j in 1:n_gs) {
    log_p[n,j] = log_sum_exp(log_p_raw[n,,j]);
    p[n,j] = exp(log_p[n,j]);
    // model site abundance
    for(s in 1:S) {
    lambda[s,n,j] = exp(log_lambda_psi[s,n] + log_p[n,j] + log_activ + log(eps_ngs[s,n,j]) + log(pbar[s,n])) .* survey_area[n];
    }
  }

// Royle-Nichols implementation in STAN (looping over possible discrete values of N)
// https://discourse.mc-stan.org/t/royle-and-nichols/14150
// https://discourse.mc-stan.org/t/identifiability-across-levels-in-occupancy-model/5340/2
    for(s in 1:S) {
if (n_survey[n] > 0) {
  vector[n_max[n,s]] lp;
    if(any_seen[s,n] == 0){ // not seen
      lp[1] = log_sum_exp(neg_binomial_2_log_lpmf(0 | log_lambda_psi[s,n], odRN), neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]));
    } else {
      lp[1] = neg_binomial_2_log_lpmf(1 | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | r[start_idx[n]:end_idx[n]]);
    }
    for (k in 2:n_max[n,s]){
      lp[k] = neg_binomial_2_log_lpmf(k | log_lambda_psi[s,n], odRN) + bernoulli_lpmf(y2[s,start_idx[n]:end_idx[n]] | 1-(1-r[start_idx[n]:end_idx[n]])^k);
    }
    lp_site[s,n] = log_sum_exp(lp);
    } else {
    lp_site[s, n] = 0;
    }
  }
    }
}

model {

  for(s in 1:S) {
    beta_psi[s] ~ normal(0, 3); // prior for intercept in poisson model
    bioregion_sd[s] ~ normal(0,2);
    bioregion_raw[s,] ~ normal(0,1);
    to_vector(eps_raw[s,,]) ~ std_normal();
    grp_sd[s] ~ normal(0, 1);
    zeta[s,] ~ normal(0, 2);
    // od
    od_sd[s] ~ normal(0,1);
    od_mu[s] ~ normal(0,1);
    od_raw[s,] ~ normal(0,1);
    // evc re
    evc_raw[s,] ~ normal(0,1);
  }
  evc_sd ~ normal(0,1);
  odRNmu ~ normal(0,1);
  beta_trans_det ~ normal(0, 2); // beta for transect detection
  beta_det ~ normal(0, 4); // prior for sigma
  activ ~ beta(bshape, bscale);  //informative prior
  log_theta ~ normal(0,2); // prior for theta

  for(n in 1:n_site) {
  for(j in 1:n_gs) {
  for(s in 1:S) {
  target += neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s, site_bioreg[n]]);
        }
  y[n,,j] ~ multinomial_logit(to_vector(log_p_raw[n,,j]));
  }
  for(s in 1:S) {
  target += lp_site[s, n];
  }
}
}

generated quantities {
  array[S, n_site,n_gs] real n_obs_pred;
  array[S, n_site, n_gs] real n_obs_true;
  array[S, n_site] real N_site;
  array[S, n_site] real N_site_pred;
  array[n_site, max_int_dist+1] real DetCurve;
  array[n_site, n_gs] real log_lik1;
  array[S, n_site, n_gs] real log_lik2;
  array[n_site, n_gs] real log_lik2_site;
  array[S, n_site] real log_lik2_species;
  array[n_site] real log_lik;
  array[n_site] real log_lik_det;
  //array[S, n_site] real Site_lambda;
  real av_gs[S];
  array[S] simplex[n_gs] eps_gs_ave;
  array[S, npc] real pred;
  //array[S, npc] real pred_trunc;
  array[S, np_reg] real Nhat_reg;
  // array[np_reg] real Nhat_reg_design;
  real Nhat[S];
  //real Nhat_trunc[S];
  real Nhat_sum;
  int trunc_counter[S];
  trunc_counter[S] = 0;

for(s in 1:S) {
  eps_gs_ave[s] = exp(zeta[s])/sum(exp(zeta[s]));
  av_gs[s] = sum(gs .* eps_gs_ave[s]); //average group size
}


for(n in 1:n_site) {
  for(j in 1:n_gs) {
  log_lik1[n,j] = multinomial_logit_lpmf(y[n,,j] | to_vector(log_p_raw[n,,j])); //for loo
  for(s in 1:S) {
  log_lik2[s,n,j] = neg_binomial_2_lpmf(n_obs[n,j,s] | lambda[s,n,j], od[s, site_bioreg[n]]); //for loo
  n_obs_true[s, n, j] = gs[j] * (neg_binomial_2_log_rng(log_lambda_psi[s,n] + log(eps_ngs[s,n,j]), od[s, site_bioreg[n]]));
  n_obs_pred[s, n,j] = gs[j] * neg_binomial_2_rng(lambda[s,n,j], od[s, site_bioreg[n]]);
  }
  log_lik2_site[n, j] = log_sum_exp(log_lik2[,n,j]);
    }
    // get loglik on a site level
    log_lik_det[n] = log_sum_exp(log_lik1[n,]);
    log_lik[n] = log_sum_exp(log_sum_exp(log_lik_det[n],
    log_sum_exp(log_lik2_site[n,])), log_sum_exp(lp_site[,n]));
      for(s in 1:S) {
    //Site_lambda[s,n] = exp(log_lambda_psi[s,n]);
    log_lik2_species[s, n] = log_sum_exp(log_sum_exp(log_lik2[s,n,]), lp_site[s,n]);
    N_site[s,n] = sum(n_obs_true[s,n,]);
    N_site_pred[s,n] = sum(n_obs_pred[s,n,]);
      }

  // loop over distance bins
  if(keyfun == 0) {
  for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
    DetCurve[n, j+1] = exp(-(j+0.5)^2/(2*sigma[n]^2)); // half normal
    }
  } else if(keyfun == 1) {
    for(j in 0:max_int_dist) { // get DS predictions for distance 0 to max bin distance
    DetCurve[n, j+1] =  1 - exp(-((j+0.5)/sigma[n])^(-theta)); //hazard rate
    }
  }
}

for(s in 1:S) {
for(i in 1:np_reg) Nhat_reg[s,i] = 0;

for(i in 1:npc) {
  pred[s,i] = neg_binomial_2_log_rng(X_pred_psi[i,] * beta_psi[s] + eps_bioregion[s, pred_bioreg[i]] + eps_evc[s, pred_evc[i]], od[s, pred_bioreg[i]]) * prop_pred[i] * av_gs[s]; //offset
  if(pred[s,i] > max(N_site[s,])) {
    //pred_trunc[s,i] = max(N_site[s,]);
    trunc_counter[s] += 1;
  }
  Nhat_reg[s,pred_reg[i]] += pred[s,i];
}
Nhat[s] = sum(pred[s,]);
//Nhat_trunc[s] = sum(pred_trunc[s,]);
}
Nhat_sum = sum(Nhat);
}

Show code

model_negbin_ms <- cmdstan_model(here::here("stan", "count_det_nondet_negbin_ms.stan"))
model_negbin_ms2 <- cmdstan_model(here::here("stan", "count_det_nondet_negbin_ms_2.stan"))
model_negbin_ms_evc <- cmdstan_model(here::here("stan", "count_det_nondet_negbin_ms_evc.stan"))

Fit models

Detection Models

Using the distance package we compare several detection functions (hazard and half-normal) with the inclusion of a parameter (herbaceous understorey cover), which may impact detection rates. The top model (as chosen by AIC) from this selection will be included in our full STAN model.

Show code

filter_behaviour <- TRUE
snapshot_interval <- 2
  # Camera information
  theta <- 40 * pi / 180 # camera angle in radians

  ### Get camera information for the camera sites
  ### For cameras that have problems you need to use the date from when the problem started
  cams_curated2 <- dplyr::tbl(con, dbplyr::in_schema("camtrap", "curated_camtrap_operation")) %>%
    dplyr::filter(ProjectShortName %in% !!project_short_name) %>% # Only retrieve cam records for hog deer 2023
    dplyr::collect() %>%
    dplyr::mutate(DateTimeDeploy = as.POSIXct(DateTimeDeploy),
                  DateTimeRetrieve = as.POSIXct(DateTimeRetrieve),
                  Tk = as.numeric(DateTimeRetrieve - DateTimeDeploy, units = "secs"), #seconds
                  Tk_prob = dplyr::coalesce(as.numeric(as.POSIXct(Problem1_to,
                                                                  format = "%Y-%m-%d %H:%M:%OS") -
                                                         as.POSIXct(Problem1_from, format = "%Y-%m-%d %H:%M:%OS"),
                                                       units = "secs"), 0),
                  Tk_adj = Tk-Tk_prob,
                  Tkt = Tk_adj / snapshot_interval, # snapshot moments: every second second
                  Effort = (Tk_adj*theta)/(snapshot_interval * pi * 2)) %>%
    dplyr::arrange(SiteID)


  ### Get the deer records (hog deer)
  camtrap_records_deer <- dplyr::tbl(con, dbplyr::in_schema(schema = "camtrap", table = "curated_camtrap_records")) %>%
    dplyr::filter(scientific_name %in% deer_species_all & ProjectShortName %in% !!project_short_name) %>%
    dplyr::select(Species = scientific_name, SiteID, Distance = metadata_Distance, size = metadata_Multiples, Date, Time, Behaviour = metadata_Behaviour) %>%
    dplyr::collect() %>%
    dplyr::mutate(Distance = dplyr::case_when(Distance == "NA" | is.na(Distance) ~ "999",
                                              TRUE ~ Distance)) %>%
    dplyr::mutate(Time = as.POSIXct(Time, format = "%H:%M:%OS"),
                  Time_n = as.numeric(Time, units = "secs"),
                  Behaviour = na_if(x = Behaviour, y = "NA")) %>%
    dplyr::rowwise() %>%
    dplyr::mutate(DistanceMod = list(stringr::str_split(Distance, pattern = "_&_")[[1]])) %>%
    dplyr::mutate(Distance = DistanceMod[which.min(as.numeric(stringr::str_extract(DistanceMod, pattern = "[0-9]+")))]) %>%
    dplyr::mutate(Distance = dplyr::case_when(Distance == "0-2.5" ~ "0 - 2.5",
                                              Distance == "999" ~ NA_character_,
                                              TRUE ~ Distance)) %>%
    dplyr::ungroup() %>%
    dplyr::filter(Time_n %% snapshot_interval == 0) %>% #& #snapshot moment interval of 2s
    {if(filter_behaviour) dplyr::filter(., is.na(.data[["Behaviour"]])) else .} %>% # filter out behaviors such as camera or marker interaction
    dplyr::group_by(SiteID, Time_n, Species) %>%
    dplyr::slice(1) %>% # if two photos occur in one second take only one (snapshot moment = 2)
    dplyr::ungroup()

  ### Tidy format for the records
  dcount<- camtrap_records_deer %>%
    dplyr::select(Species, SiteID, Distance, size, Date, Time, Behaviour) %>%
    dplyr::full_join(cams_curated2 %>%
                       dplyr::select(SiteID, Tkt, Effort), by="SiteID") %>%
    dplyr::mutate(object=1:nrow(.)) %>%
    dplyr::mutate(size = if_else(is.na(size),0L, size)) %>%
    dplyr::arrange(SiteID)

  ### Format data in a time format for availability
  summarised_count <- dcount
  summarised_count$Distance <- factor(summarised_count$Distance, levels = c("0 - 2.5", "2.5 - 5", "5 - 7.5", "7.5 - 10", "10+"))
  
  summarised_count <- summarised_count %>%
    mutate(distance = case_when(Distance == "0 - 2.5" ~ 1.25,
                                Distance == "2.5 - 5" ~ 3.75,
                                Distance == "5 - 7.5" ~ 6.25,
                                Distance == "7.5 - 10" ~ 8.75,
                                Distance == "10+" ~ 11.25)) %>%
    filter(size == 1)
  
  # site variables 
    site_vars <- cams_curated2 %>%
    left_join(dplyr::tbl(con, dbplyr::in_schema("deervic", "curated_site_data")) %>%
    dplyr::filter(SiteID %in% !!c(cams_curated2$SiteID, "47191")) %>%
    dplyr::collect() %>%
    dplyr::mutate(HerbaceousUnderstoryCover = scale(NNWHUCover + ENWHUCover),
                  WoodyUnderstoryCover = scale(NWUCover + EWUCover),
                  Trail = factor(Trail),
                  SiteID = dplyr::case_when(SiteID == "47191" & CameraID == "HO04101053" ~ "47191A",
                                            SiteID == "47191" & CameraID != "HO04101053" ~ "47191B",
                                            TRUE ~ SiteID))) %>% # native + exotic herbaceous cover
    dplyr::arrange(SiteID) %>%
    as.data.frame() %>%
      select(-Time, -Tkt, -Date, -Effort) 
    
    ds_data <- summarised_count %>% left_join(site_vars) %>%
      filter(!is.na(HerbaceousUnderstoryCover))
  
  mybreaks <- c(0,2.5,5,7.5,10,12.5)
  trunc.list <- list(left=0, right=12.5)
  conversion <- convert_units("metre", "metre", "metre")
  
  hn0 <- ds(ds_data, transect = "point", key="hn", adjustment = NULL,
          convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
  
  hr0 <- ds(ds_data, transect = "point", key="hr", adjustment = NULL,
          convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
  
  hn1 <- ds(ds_data, transect = "point", key="hn", adjustment = NULL, 
            formula = ~HerbaceousUnderstoryCover,
          convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
  
  hr1 <- ds(ds_data, transect = "point", key="hr", adjustment = NULL, 
            formula = ~HerbaceousUnderstoryCover,
          convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list)
  
  kableExtra::kbl(summarize_ds_models(hn0, hr0, hn1, hr1, output = "latex") %>% 
                    mutate(Model = stringr::str_extract(Model, pattern = '(?<=\\{)[^\\}]+')), 
                  digits=3,
             caption="Model selection table for deer detection", format = "html") %>%
    kable_styling(bootstrap_options = "striped")

Table 1: Model selection table for deer detection

	Model	Key function	Formula	\(\chi^2\) \(p\)-value	\(\hat{P_a}\)	se(\(\hat{P_a}\))	\(\Delta\)AIC
4	hr1	Hazard-rate	~HerbaceousUnderstoryCover	0	0.230	0.008	0.000
2	hr0	Hazard-rate	~1	0	0.228	0.008	3.709
1	hn0	Half-normal	~1	0	0.336	0.004	186.819
3	hn1	Half-normal	~HerbaceousUnderstoryCover	0	0.336	0.004	187.766

Model Options

Below we list the models we are implementing, all have a hazard detection function, but we compare four different abundance formulas as well as the inclusion of the EVC random effect.

Show code

#Models to run 
models_to_run <- data.frame(formula = c("ab_formula_1", 
                                        "ab_formula_2",
                                        "ab_formula_3",
                                        "ab_formula_4",
                                        "ab_formula_3", 
                                        "ab_formula_4",
                                        "ab_formula_1",
                                        "ab_formula_2"), 
                            keyfun = c(1,1,1,1,1,1,1,1), 
                            evc = c(F,F,F,F,T,T,T,T))

Full Models

We can fit models using cmdstanr. Here we fit models using a poisson distribution. The models we compare are:

• Model with several key abundance covariates and half-normal detection function
• Model with key abundance covariates and hazard-rate detection function
• Null-model (bioregion only) with half-normal detection function
• Model with key abundance covariates, bioregion random-effect and an EVC group random effect.

Show code

model_fits <- list()

for(i in 1:nrow(models_to_run)) {

  if(models_to_run$evc[i]) {
    model_to_fit <- model_negbin_ms_evc
  } else {
    model_to_fit <- model_negbin_ms
  }
  
  form_to_use <- get(models_to_run$formula[i])
  
  params_to_use <- c(TRUE, labels(terms(ab_formula_1)) %in% labels(terms(form_to_use)))
  
  data_to_use <- multispecies_data
  
  data_to_use$X_psi <- as.matrix(data_to_use$X_psi[,params_to_use])
  data_to_use$X_pred_psi <- as.matrix(data_to_use$X_pred_psi[,params_to_use])
  data_to_use$m_psi <- sum(params_to_use)
  
  data_to_use$keyfun <- models_to_run$keyfun[i]

  model_fits[[i]] <- model_to_fit$sample(data = data_to_use,
                                         chains = nc,
                                 parallel_chains = nc, 
                                 init = 0.1, 
                                 max_treedepth = 10, 
                                 refresh = 25, 
                                 show_messages = TRUE,
                                 show_exceptions = FALSE,
                                 save_warmup = FALSE,
                                 iter_sampling = ni, 
                                 iter_warmup = nw)

  model_fits[[i]]$save_object(paste0("outputs/models/fit_multispecies_", 
                                     i,".rds"))

  
}

Show code

prefix <- "outputs/models/fit_multispecies_"
model_fits <- list()
for(i in 1:nrow(models_to_run)) {
  file_to_read <- paste0(prefix, i, ".rds")
  
  model_fits[[i]] <- readRDS(file_to_read)
}

Model Evaluation

Our strategy for model evaluation is to determine the best performing model out of a limited set of options with different fixed-effects, random-effects and detection functions.

LOO (Leave-One-Out Cross-Validation).

We use leave-one-out cross-validation (LOO-CV); a Bayesian model evaluation process (Vehtari, Gelman, and Gabry 2017; Vehtari et al. 2020) to compare the models. Better performing models according to loo::loo() will have higher elpd values.

Show code

loos <- list()
for(i in 1:length(model_fits)) {
  loos[[i]] <- model_fits[[i]]$loo(variables = "log_lik2_species", cores = 8)
}

names(loos) <- paste("model", 1:length(model_fits))

loo_compare_table <- loo::loo_compare(x = loos)

# Plot loo table
  gt(loo_compare_table %>% 
       as.data.frame() %>% 
              tibble::rownames_to_column()) %>% 
    fmt_number(decimals = 2) %>%
    tab_style(locations = cells_row_groups(), style = cell_text(weight = "bold"))

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic
model 7	0.00	0.00	2,120.49	7.69	1.52	0.08	−4,240.97	15.38
model 8	−0.42	0.35	2,120.07	7.69	1.46	0.08	−4,240.13	15.37
model 5	−0.56	0.38	2,119.92	7.67	1.43	0.08	−4,239.85	15.34
model 6	−2.15	0.76	2,118.34	7.67	1.28	0.07	−4,236.68	15.33
model 1	−2.73	0.74	2,117.75	7.70	1.27	0.07	−4,235.51	15.40
model 2	−3.14	0.82	2,117.35	7.70	1.21	0.07	−4,234.70	15.39
model 3	−3.39	0.84	2,117.10	7.69	1.18	0.07	−4,234.20	15.38
model 4	−5.70	1.05	2,114.79	7.70	1.05	0.06	−4,229.57	15.39

Given that models are similar in there predictive performance at sampled sites we also evaluate the performance of predictions across unsampled sites. Below we see that more complex models (usually having a slighly better elpd), have poorer performance to unsampled areas and lead to unrealistic predictions.

Show code

Nhat_CV <- list()
for(i in 1:length(model_fits)) {
  Nhat_CV[[i]] <- model_fits[[i]]$summary("Nhat", CV = ~ sd(.x)/mean(.x)) %>%
    mutate(Model = i, Species = deer_species_all)
}

Nhat_CV_c <- bind_rows(Nhat_CV)

Nhat_CV_c %>%
  select(Model, Species, CV) %>%
  kbl("html", escape=FALSE, caption = "Coefficients of Variation for Total Abundance estimates for the models assessed in the analysis") %>%
  kable_styling() %>%
  kableExtra::column_spec(3, background = ifelse(Nhat_CV_c$CV > 0.5, "red", ifelse(Nhat_CV_c$CV > 0.2, "orange", "green")), color = "white") %>%
 collapse_rows(1) 

Table 2: Coefficients of Variation for Total Abundance estimates for the models assessed in the analysis

Model	Species	CV
1	Cervus unicolor	0.1530477
	Dama dama	0.4989131
	Cervus elaphus	0.6918822
	Axis porcinus	0.4821552
2	Cervus unicolor	0.1674774
	Dama dama	0.3884288
	Cervus elaphus	0.4659085
	Axis porcinus	0.3418980
3	Cervus unicolor	0.1391707
	Dama dama	0.3201355
	Cervus elaphus	0.3513872
	Axis porcinus	0.3565110
4	Cervus unicolor	0.1368803
	Dama dama	0.1942815
	Cervus elaphus	0.3120815
	Axis porcinus	0.2593997
5	Cervus unicolor	0.1829863
	Dama dama	0.3010461
	Cervus elaphus	1.3648098
	Axis porcinus	0.5658459
6	Cervus unicolor	0.2034979
	Dama dama	0.2388721
	Cervus elaphus	0.5028944
	Axis porcinus	0.7388782
7	Cervus unicolor	0.1985761
	Dama dama	0.3640237
	Cervus elaphus	4.4050792
	Axis porcinus	1.0583120
8	Cervus unicolor	0.1899283
	Dama dama	0.3565456
	Cervus elaphus	1.2353806
	Axis porcinus	0.5980698

Selecting the ‘top’ model

Based on the findings of loo and comparisons of other checks (posterior predictive checks, variation in posterior draws, plausibility of predictions), we see that all model perform similarly well using loo. However, several models produce predictions with lower levels of variation in the posterior predictions. Based on this, we use model 3 hereafter for analysis. This model has fixed effects of scale(sqrt(BareSoil)), scale(sqrt(Nitrogen)), scale(log(PastureDistance)), scale(BIO15) and scale(sqrt(ForestEdge)). It also has a random effect of bioregion for each species and an over-dispersion parameter in the camera counts based on bioregion. This model does not have as many fixed effects as models 1 and 2, and also does not have a second random effect of EVC alongside bioregion (model 6). Unfortunately, the additional random effect of EVC produces much more variable and over-dispersed results for Red Deer, so we opt to use model 3 instead.

Show code

# assign index for top model
top <- 3

Posterior predictive checks

Posterior predictive checks allow us to compare the observed data to the model-generated data. For each species we undertake posterior predictive checks for summary statistics relating to the number of deer seen on the cameras at each site. Ideally a well-fit model is able to make predictions that match the observed data. Here counts are the number of snapshot moments with deer. The summary statistics we use for the posterior predictive checks are:

1. Maximum counts of deer seen at a site on a camera
2. Mean number of deer seen at a site on a camera
3. Standard deviation of the counts of deer on the cameras
   
4. Average (mean) counts at sites in a scatter plot
5. Proportion of sites with zero counts (camera or transect)
   
6. Observed vs expected proportion of counts

Sambar Deer

Based on the plot below, we see that the PPCs for Sambar perform sufficiently well.

Show code

q95 <- function(x) quantile(x, 0.95, na.rm = T)
# q25 <- function(x) quantile(x, 0.25, na.rm = T)
q75 <- function(x) quantile(x, 0.75, na.rm = T)
sd_90 <- function (x, na.rm = FALSE) {
  quants <- quantile(x, c(0.05, 0.95), na.rm = T)
  x <- x[x < quants[2] & x > quants[1]]
  sqrt(var(if (is.vector(x) || is.factor(x)) x else as.double(x),
    na.rm = na.rm))
}

bayes_theme <- c(delwp_cols_seq$Teal[6:8], delwp_cols_shades$Navy[1:3])

bayesplot::color_scheme_set(bayes_theme)

funs <- c(max, mean, sd, ppc_scatter_avg, prop_zero, ppc_dens_overlay)
titles <- c("Max", "Mean", "SD", "Average Site Counts", "Proportion Zeros", "Density of Counts")

ppc_sambar <- list()

# funs <- c(prop_zero, mean, q90, sd)

for(i in 1:length(funs)) {
ppc_sambar[[i]] <- posterior_checks_multispecies(model = model_fits[[top]], 
                 model_data = multispecies_data, species_index = 1,
                 stat = funs[[i]], integrated = T,
                 title = titles[i])
}

cowplot::plot_grid(plotlist = ppc_sambar, labels = "AUTO", ncol = 2)

Figure 1: Posterior predictive checks for Sambar Deer

Fallow Deer

The PPC’s below show that fallow has more dispersion in camera counts than Sambar. While the observed counts fall within the intervals of our model, there is higher uncertainty. Some of these findings appear to be driven by a site with very high counts (over 2000 snapshot moments of deer individuals), this is more than twice the second highest.

Show code

ppc_fallow <- list()

for(i in 1:length(funs)) {
ppc_fallow[[i]] <- posterior_checks_multispecies(model = model_fits[[top]], 
                 model_data = multispecies_data, species_index = 2,
                 stat = funs[[i]], integrated = T,
                 title = titles[i])
}

cowplot::plot_grid(plotlist = ppc_fallow, labels = "AUTO", ncol = 2)

Figure 2: Posterior predictive checks for Fallow Deer

Red Deer

Red Deer have PPC’s that show a high congruence between predicted and expected counts.

Show code

ppc_red <- list()

for(i in 1:length(funs)) {
ppc_red[[i]] <- posterior_checks_multispecies(model = model_fits[[top]], 
                 model_data = multispecies_data, species_index = 3,
                 stat = funs[[i]], integrated = T, 
                 title = titles[i])
}

cowplot::plot_grid(plotlist = ppc_red, labels = "AUTO", ncol = 2)

Figure 3: Posterior predictive checks for Red Deer

Hog Deer

Hog Deer have PPC’s that show our model performs sufficiently well in explaining the patterns of observed counts.

Show code

ppc_hog <- list()

for(i in 1:length(funs)) {
ppc_hog[[i]] <- posterior_checks_multispecies(model = model_fits[[top]], 
                 model_data = multispecies_data, species_index = 4,
                 stat = funs[[i]], integrated = T,
                 title = titles[i])
}

cowplot::plot_grid(plotlist = ppc_hog, labels = "AUTO", ncol = 2)

Figure 4: Posterior predictive checks for Hog Deer

Model Convergence

We observed no divergences or any STAN sampling warnings/issues for our models. To visualise the convergence of the model we can observe the mixing of chains for key parameters below:

• beta_det: Parameters for the distance sampling model (intercept and herbaceous understorey)
  
• beta_psi: Parameters for the fixed effect abundance model
  
• bioregion_sd: Parameter for the bioregion random effect standard deviation
  
• odRN: Parameter for the over-dispersion in the Royle-Nichols negative binomial model
  
• od_sd and od_sd: Parameters for the mean and standard deviation of over-dispersion link with bioregion for the camera counts.

Show code

bayesplot::color_scheme_set(unname(delwp_cols[1:6]))
convergence_params <- model_fits[[top]]$draws(c("beta_det", 
                                                "beta_psi",
                                                "bioregion_sd", 
                                                "beta_trans_det", 
                                                "odRN", 
                                                "od_mu", 
                                                "od_sd"))

mcmc_trace(convergence_params, facet_args = list(ncol = 4)) + 
  theme(panel.spacing = unit(0.5, "lines"))

Figure 5: Chain mixing of several key parameters relating to detection, abundance and dispersion

Model Predictions

Within the STAN model we generate predictions for sampled and un-sampled locations. This provides us with site-level abundance estimates as well as estimates across all (un-sampled) public forest.

Site-based Predictions

Site-based Prediction Map

Visualisation of point-estimates for the various deer species surveyed for in this study.

Show code

vic_regions <- vicmap_query("open-data-platform:delwp_region") %>%
  collect() %>%
  st_transform(3111) %>%
  st_simplify(dTolerance = 500)

delwp_pal <- colorRampPalette(c("#B9C600",
                                delwp_cols[["Teal"]], 
                                delwp_cols[["Navy"]]))


site_preds <- function(model, cams_curated, species_index) {
  
rn_dens <- model$summary("N_site", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), posterior::default_summary_measures(), posterior::default_convergence_measures())
  
which_sp <- which(stringr::str_detect(string = rn_dens$variable,
                                         pattern = paste0("N_site\\[",
                                                          species_index)))
  

  density_at_sites_rn <- cbind(cams_curated, rn_dens[which_sp, ])
  
  return(density_at_sites_rn)
}

site_density_plot <- function(densities, regions, species) {
  densities$density <- cut(densities$mean,
                                   breaks = c(0, 0.25, 0.5, 1, 3, 5, 10, max(densities$mean)),
                                   labels = c("< 0.25",
                                                    "0.25 - 0.5", 
                                                    "0.5 - 1", 
                                                    "1 - 3", 
                                                    "3 - 5",
                                                    "5 - 10", 
                                                    "10 +"), include.lowest = T, right = T)
  
  if(species == "Hog") {
    regions <- regions %>% filter(delwp_region == "GIPPSLAND")
    densities <- densities %>% st_filter(regions %>% st_transform(4283))
  }
  
   lvls <- length(unique(densities$density, na.rm = T))
  
plot <- ggplot2::ggplot(data = densities) +
  ggplot2::geom_sf(data = regions, alpha = 0.75, fill = "grey80") +
  ggplot2::geom_sf(aes(fill = density, alpha = mean), shape = 21, size = 3) +
  ggplot2::scale_fill_manual(values = delwp_pal(lvls), 
                             na.value = "transparent", na.translate = F,
                             name = "", guide = guide_legend(override.aes = list(size = 6))) +
  scale_alpha_continuous(range = c(0.5,1), guide = "none") +
  labs(title = paste0('Abundance of ',species ,' Deer'), fill = bquote('Deer per'~km^2)) +
  theme_bw() +
  theme(legend.text = element_text(size = 18), legend.key.size = unit(1, "cm"),
        title = element_text(size = 20))
  
return(plot)
  
}

sambar_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 1)
fallow_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 2)
red_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 3)
hog_preds <- site_preds(model_fits[[top]], cams_curated = cams_curated, species_index = 4)

site_deer_predictions <- list(Sambar = sambar_preds, 
                         Fallow = fallow_preds, 
                         Red = red_preds, 
                         Hog = hog_preds)

saveRDS(site_deer_predictions, "outputs/site_deer_predictions.rds")

cowplot::plot_grid(site_density_plot(sambar_preds, vic_regions, species = "Sambar"),
site_density_plot(fallow_preds, vic_regions, species = "Fallow"),
site_density_plot(red_preds, vic_regions, species = "Red"), 
site_density_plot(hog_preds, vic_regions, species = "Hog"), ncol = 1)

Figure 6: Site-level density estimates for all sites sampled as part of statewide and hog deer surveys. Point-density estimates are from trimmed means of the posterior draws (trimming of the top and bottom 1% of draws)

Site-density Summaries

As a sanity-check we compare the average modelled densities at sites with (i) no evidence of deer, (ii) evidence of deer present on transects, and (iii) evidence of deer present on cameras. We expect that average densities are generally higher at sites that detected some form of deer than sites that did not detect any sign of deer. Additionally, we would also expect average densities to be generally higher at sites that had detections on cameras than those with only detections from transects. The table below shows these expectations to be correct.

Show code

density_summary_table <- function(preds, model_data, species, species_index) {
  
  cam_seen <- as.integer(as.logical(rowSums(model_data$n_obs[,,species_index])))
  
  preds_sum <- preds %>% 
    st_drop_geometry() %>%
    mutate(`Species` = species,
           `CameraDetection` = cam_seen, 
           `AnyDetection` = model_data$any_seen[species_index, ], 
           `Detection` = case_when(CameraDetection == 0 & AnyDetection == 0 ~ "Not seen", 
                                   CameraDetection == 0 & AnyDetection == 1 ~ "Only detected on transects", 
                                   CameraDetection == 1 & AnyDetection == 1 ~ "Seen on cameras")) %>%
    group_by(`Species`, `Detection`) %>%
    summarise(`Number of Sites` = n(),
              `Average Density` = mean(trimmed_mean)) %>%
    ungroup()
  
  return(preds_sum)
}

density_summary <- bind_rows(
  density_summary_table(sambar_preds, multispecies_data, 
                        species = "Sambar", species_index = 1),
  density_summary_table(fallow_preds, multispecies_data, 
                        species = "Fallow", species_index = 2),
  density_summary_table(red_preds, multispecies_data, 
                        species = "Red", species_index = 3),
  density_summary_table(hog_preds, multispecies_data, 
                        species = "Hog", species_index = 4))

density_summary %>%
  kbl(format = "html", digits = 2, caption = "Average (98% trimmed mean) density estimates at the various sites groups of sites") %>%
  kable_styling("striped") %>%
  collapse_rows(1)

Table 3: Average (98% trimmed mean) density estimates at the various sites groups of sites

Species	Detection	Number of Sites	Average Density
Sambar	Not seen	180	0.81
	Only detected on transects	33	2.31
	Seen on cameras	104	3.62
Fallow	Not seen	253	0.28
	Only detected on transects	34	0.31
	Seen on cameras	30	1.02
Red	Not seen	304	0.03
	Only detected on transects	1	1.22
	Seen on cameras	12	1.38
Hog	Not seen	295	0.23
	Seen on cameras	22	2.62

Regional and Statewide Abundance

Within our model we calculate abundance/density for deer in each of the 118947 km² of public land. Based on these spatial predictions we can estimate abundance at a regional level (6 DEECA regions) and across the whole state.

Statewide Maps

Using the model predictions (“pred”) for all suitable public land we generate a raster (1km² resolution). We save the average spatial estimates under outputs/rasters (tif files) and also provide binned plots (png files) in outputs/plots.

Show code

pred_raster_full <- terra::rast(prediction_raster)

pred_raster <- terra::app(pred_raster_full[[stringr::str_subset(
  stringr::str_remove_all(labels(terms(ab_formula_1)),
                          "scale[(]|[)]|log[(]|sqrt[(]"), 
  pattern = "[*]", negate = T)]], mean)

mean_raster <- list()
occ_raster <- list()
mean_raster_discrete <- list()

# gp_preds_draws_all <- model_fits[[top]]$draws("pred", format = "matrix")
gp_preds_draws_all <- model_fits[[top]]$draws("pred", format = "matrix")
# gp_preds_summary_all <- model_fits[[top]]$summary("pred", prob_occ = ~ mean(. > 0))


for(i in 1:length(deer_species_all)) {

which_sp <- which(stringr::str_detect(string = colnames(gp_preds_draws_all),
                                         pattern = paste0("pred\\[", i)))

gp_preds_draws_sp <- gp_preds_draws_all[,which_sp]

terra::values(pred_raster)[!is.na(terra::values(pred_raster))] <- apply(gp_preds_draws_sp, 
                                                                        2, 
                                                                        mean, 
                                                                        na.rm = T, 
                                                                        trim = 0.001) #filter out highest draw

mean_raster[[i]] <- pred_raster
mean_raster_discrete[[i]] <- mean_raster[[i]]
max_pred <- max(values(mean_raster[[i]]), na.rm = T)

values(mean_raster_discrete[[i]]) <- cut(values(mean_raster_discrete[[i]]) , 
                                         breaks = c(0,0.25,0.5,1,3,5,10,max_pred), 
                                         include.lowest = T, right = T,
                                         labels = c("< 0.25",
                                                    "0.25 - 0.5", 
                                                    "0.5 - 1", 
                                                    "1 - 3", 
                                                    "3 - 5",
                                                    "5 - 10", 
                                                    "10 +"))

occ_raster[[i]] <- pred_raster
terra::values(occ_raster[[i]])[!is.na(terra::values(occ_raster[[i]]))] <- apply(gp_preds_draws_sp, 2, function(x) sum(x > 0)/length(x))

}

# combine mean rasters together
combined_raster <- rast(mean_raster)
names(combined_raster) <- c("Average Sambar Deer Density (per km2)", 
                            "Average Fallow Deer Density (per km2)", 
                            "Average Red Deer Density (per km2)", 
                            "Average Hog Deer Density (per km2)")

combined_occupancy_raster <- rast(occ_raster)
names(combined_occupancy_raster) <- c("Average Sambar Deer Occupancy Probability", 
                            "Average Fallow Deer Occupancy Probability", 
                            "Average Red Deer Occupancy Probability", 
                            "Average Hog Deer Occupancy Probability")

writeRaster(combined_occupancy_raster, "outputs/rasters/combined_occupancy_raster.tif", overwrite = T)
writeRaster(combined_raster, "outputs/rasters/combined_deer_average_density.tif", overwrite = T)

Show code

# reg <- VicmapR::vicmap_query("open-data-platform:delwp_region") %>%
#     VicmapR::collect() %>%
#   sf::st_transform(3111) %>%
#   sf::st_simplify(dTolerance = 250) 
state <- VicmapR::vicmap_query("open-data-platform:vmlite_victoria_polygon_su5") %>%
  filter(state != "NSW" & state != "SA" & feature_type_code != "sea") %>%
  VicmapR::collect() %>%
  sf::st_transform(3111)

gippsland <- vic_regions %>% filter(delwp_region == "GIPPSLAND")

plot_abundance <- function(raster, state, species, crop = NULL) {
  
  if(!is.null(crop)) {
    raster <- terra::crop(raster, vect(crop), mask = T)
    state <- crop
  }
  
  lvls <- length(unique(values(raster, na.rm = T)))
  
  plot <- ggplot2::ggplot() +
    ggplot2::geom_sf(data = state, alpha = 1, linewidth = 0.5, fill = "grey90") + 
    tidyterra::geom_spatraster(data = raster, na.rm = T) + 
    # tidyterra::scale_fill_terrain_d(na.translate = FALSE) + 
    ggplot2::scale_fill_manual(values = delwp_pal(lvls), na.value = "transparent", na.translate = F) +
    ggplot2::labs(fill = bquote('Deer per'~km^2), title = paste0("Average Density of ", species, " on Victorian Public Land")) + 
    ggspatial::annotation_scale() +
    delwp_theme()
  
  return(plot)
  
} 

sambar_abundance_plot <- plot_abundance(mean_raster_discrete[[1]], 
                                        state = state, 
                                        species = "Sambar Deer")
fallow_abundance_plot <- plot_abundance(mean_raster_discrete[[2]], 
                                        state = state, 
                                        species = "Fallow Deer")
red_abundance_plot <- plot_abundance(mean_raster_discrete[[3]], 
                                     state = state, 
                                     species = "Red Deer")
hog_abundance_plot <- plot_abundance(mean_raster_discrete[[4]], 
                                     state = state, crop = gippsland,
                                     species = "Hog Deer")

ggsave(plot = sambar_abundance_plot, filename = "outputs/plots/sambar_abundance_plot.png", 
       width = 2400, height = 1600, units = "px")
ggsave(plot = fallow_abundance_plot, filename = "outputs/plots/fallow_abundance_plot.png", 
       width = 2400, height = 1600, units = "px")
ggsave(plot = red_abundance_plot, filename = "outputs/plots/red_abundance_plot.png", 
       width = 2400, height = 1600, units = "px")
ggsave(plot = hog_abundance_plot, filename = "outputs/plots/hog_abundance_plot.png", 
       width = 2400, height = 1600, units = "px")

cowplot::plot_grid(sambar_abundance_plot, fallow_abundance_plot, red_abundance_plot, hog_abundance_plot, ncol = 1)

Figure 7: Abundance of Sambar, Fallow, Red and Hog Deer across Victoria, dark-grey area reflects area not included in predictions (i.e. not public land)

Regional Abundance Estimates

For each DEECA region we provide species-level estimates of abundance with 90 % confidence interval. We also calculate the model-based average density within each region based on the abundance and the total area of public land within each region.

Show code

reg <- VicmapR::vicmap_query("open-data-platform:delwp_region") %>%
    VicmapR::collect() %>%
    dplyr::group_by(delwp_region) %>%
    dplyr::summarise(geometry = sf::st_combine(geometry)) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(delwp_region_fact = as.integer(factor(delwp_region))) 

abundance_table <- function(model, regions, pred_area, caption, species_index) {
  
  reg <- model$summary("Nhat_reg", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), posterior::default_summary_measures(), posterior::default_convergence_measures())
  
  which_sp <- which(stringr::str_detect(string = reg$variable,
                                         pattern = paste0("Nhat_reg\\[", species_index)))
  
  regional_abundance <- reg[which_sp,] %>%
    dplyr::bind_rows(model$summary("Nhat", prob_occ = ~ mean(. > 0), trimmed_mean = ~mean(., trim = 0.025), posterior::default_summary_measures(), posterior::default_convergence_measures())[species_index,]) %>% 
    dplyr::mutate(variable = c(regions$delwp_region, "TOTAL"), 
                  `Area km2` = round(c(pred_area, sum(pred_area))), 
                  `Average Density (per km2)` = round(median/`Area km2`, 2)) %>% 
    dplyr::transmute(Region = variable, 
                  `Mean` = round(mean), 
                  Median = round(median), 
                  SD = round(sd), 
                  MAD = round(mad),
                  `5 %` = round(q5), 
                  `95 %` = round(q95), 
                  `Area km2`,
                  `Average Density (per km2)`)
  
  kableExtra::kbl(regional_abundance, digits = 2, format = "html", caption = caption, format.args = list(big.mark = ",")) %>%
    kableExtra::kable_styling("striped") %>%
    kableExtra::column_spec(1, bold = TRUE) %>%
    kableExtra::row_spec(6, hline_after = T) %>%
    kableExtra::row_spec(7, background = "#A5A1B5", color = "black", bold = T, hline_after = T)
}

pred_a_df <- data.frame(pred_reg = multispecies_data$pred_reg, 
                        area = multispecies_data$prop_pred) %>%
  group_by(pred_reg) %>%
  summarise(area = sum(area))

abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Sambar Deer Density", species_index = 1)

Table 4: Regional estimates of Sambar Deer Density

Region	Mean	Median	SD	MAD	5 %	95 %	Area km2	Average Density (per km2)
BARWON SOUTH WEST	441	351	378	207	116	1,038	4,745	0.07
GIPPSLAND	79,111	77,855	11,987	11,512	61,310	100,252	25,211	3.09
GRAMPIANS	3,167	2,827	2,450	891	1,754	5,244	9,896	0.29
HUME	79,089	77,473	13,302	12,438	59,966	102,457	16,890	4.59
LODDON MALLEE	3,058	2,426	2,682	987	1,343	6,335	15,597	0.16
PORT PHILLIP	7,525	7,309	1,643	1,546	5,262	10,636	2,231	3.28
TOTAL	172,390	170,427	23,992	23,450	137,276	214,214	74,570	2.29

Show code

abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Fallow Deer Density", species_index = 2)

Table 4: Regional estimates of Fallow Deer Density

Region	Mean	Median	SD	MAD	5 %	95 %	Area km2	Average Density (per km2)
BARWON SOUTH WEST	19,587	15,193	16,262	8,464	6,347	48,250	4,745	3.20
GIPPSLAND	16,774	15,915	4,867	4,220	10,524	25,457	25,211	0.63
GRAMPIANS	7,909	6,435	9,453	2,007	3,949	13,774	9,896	0.65
HUME	24,926	23,767	7,045	6,168	15,735	37,570	16,890	1.41
LODDON MALLEE	4,573	3,913	2,606	1,465	2,269	9,293	15,597	0.25
PORT PHILLIP	2,296	2,162	851	677	1,265	3,792	2,231	0.97
TOTAL	76,065	71,419	24,351	16,799	48,994	116,448	74,570	0.96

Show code

abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Red Deer Density", species_index = 3)

Table 4: Regional estimates of Red Deer Density

Region	Mean	Median	SD	MAD	5 %	95 %	Area km2	Average Density (per km2)
BARWON SOUTH WEST	4,979	4,031	3,749	1,661	2,235	10,594	4,745	0.85
GIPPSLAND	1,015	931	459	374	438	1,892	25,211	0.04
GRAMPIANS	4,021	3,716	1,597	1,256	2,154	6,862	9,896	0.38
HUME	3,353	3,086	1,502	1,300	1,461	5,997	16,890	0.18
LODDON MALLEE	433	168	862	193	13	1,717	15,597	0.01
PORT PHILLIP	177	149	124	78	58	377	2,231	0.07
TOTAL	13,977	13,055	4,911	3,440	8,735	21,839	74,570	0.18

Show code

abundance_table(model_fits[[top]], regions = reg, pred_area = pred_a_df$area, caption = "Regional estimates of Hog Deer Density", species_index = 4)

Table 4: Regional estimates of Hog Deer Density

Region	Mean	Median	SD	MAD	5 %	95 %	Area km2	Average Density (per km2)
BARWON SOUTH WEST	153	76	416	82	6	461	4,745	0.02
GIPPSLAND	6,296	5,847	2,266	1,835	3,618	10,396	25,211	0.23
GRAMPIANS	122	53	450	55	4	380	9,896	0.01
HUME	99	57	122	63	3	333	16,890	0.00
LODDON MALLEE	106	24	724	30	1	316	15,597	0.00
PORT PHILLIP	1,188	1,029	660	474	479	2,339	2,231	0.46
TOTAL	7,964	7,412	2,839	2,275	4,611	13,183	74,570	0.10

Model Covariates

Our model uses covariates to inform three separate sub-models:

• A distance-sampling detection process in front of the camera
  
• An imperfect detection/observation model based on the transect methods and the camera survey
  
• An abundance model based on spatially-derived covariates, as well as a bioregion random effect

The following section explores the effect of these covariates and sub-models on detection and abundance.

Detection Processes

At a site, there are two processes that may lead us to not record deer on a camera when in reality they occupy the site (whereby the site is the home range area surrounding the camera.) Firstly, deer may be available for detection and enter the sampling area in front of the camera. However, due to a function of their distance from the camera we fail to detect this individual. Secondly, they may occupy a site (known from transect signs) but never enter the camera field of view, in this case we estimate the various detection rates from the various methods. The detection rate of the camera in this method is regarded as the spatial availability of deer for camera trap sampling.

Distance-sampling

Below we show the average detection rate for a given group of deer in front of the camera (up to 12.5m).

Show code

det_rates <- model_fits[[top]]$summary("p") %>%
  mutate(var = stringr::str_extract(variable, "(?<=\\[).+?(?=\\])")) %>%
  separate(var, into = c("site", "Group Size"), sep = ",")

av_det_rates <- det_rates %>%
  group_by(`Group Size`) %>%
  summarise(mean = mean(mean)) %>%
  ungroup() %>% 
  transmute(`Group Size`, 
            `Average Site Detection Probability` = mean)


av_det_rates %>%
  bind_rows() %>%
  arrange(`Group Size`) %>%
  kbl("html", digits = 3, caption = "Average detection rates for different group sizes") %>%
  kable_styling("striped") %>%
  collapse_rows(1)

Table 5: Average detection rates for different group sizes

Group Size	Average Site Detection Probability
1	0.306
2	0.390
3	0.451
4	0.498
5	0.537

We can generate detection function plots that show how detection rates fall-off over varying distances. Initial exploration showed that there may be some relationship between herbaceous understorey and detection probability for Sambar Deer. However, when pooled together we see that herbaceous understorey appears to have no effect on detection.

Show code

site_vars <- dplyr::tbl(con, dbplyr::in_schema("deervic", "curated_site_data")) %>%
    dplyr::filter(SiteID %in% !!c(cams_curated$SiteID, "47191")) %>%
    dplyr::collect() %>%
    dplyr::mutate(HerbaceousUnderstoryCover = NNWHUCover + ENWHUCover,
           SiteID = dplyr::case_when(SiteID == "47191" & CameraID == "HO04101053" ~ "47191A",
                              SiteID == "47191" & CameraID != "HO04101053" ~ "47191B",
                              TRUE ~ SiteID)) %>% # native + exotic herbaceous cover
    dplyr::arrange(SiteID) %>%
    as.data.frame()

  n_distance_bins <- 5
  delta <- 2.5
  midpts <- c(1.25, 3.75, 6.25, 8.75, 11.25)
  max_distance <- 12.5
  gs <- 1:5
  
  pa <- matrix(data = NA, nrow = max(gs), ncol = length(midpts))
  for(j in 1:max(gs)) {
    pa[j, ] <- sapply(midpts, function(x) {pr_closest(x-(0.5*delta),
                                                      x+(0.5*delta),
                                                      max_distance = max_distance,
                                                      n = j)})
  }


det_curve <- model_fits[[top]]$draws("DetCurve", format = "draws_matrix") %>%
  as.data.frame() %>%
  # head(250) %>%
  pivot_longer(cols = everything())

det_curve_wr <- det_curve %>%
  mutate(var = stringr::str_extract(name, "(?<=\\[).+?(?=\\])")) %>%
  separate(var, into = c("s", "Distance"), sep = ",") %>%
  mutate(Distance = as.numeric(Distance)-0.5)

det_vars_pred <- site_vars %>%
  mutate(s = as.character(1:nrow(.)),
         herbaceouslvl = cut(HerbaceousUnderstoryCover,
                             breaks = c(0, 25, 50, 75, 105), 
                             labels = c("0 - 25%", 
                                        "25 - 50%", 
                                        "50 - 75%", 
                                        "75 - 100%"),
                             include.lowest = T, right = FALSE))

det_curve_sum <- det_curve_wr %>%
  mutate(Distance = as.numeric(Distance)) %>%
  left_join(det_vars_pred) %>%
  group_by(herbaceouslvl, Distance) %>%
  summarise(median = quantile(value, 0.5),
            q5 = quantile(value, 0.05),
            q95 = quantile(value, 0.95))

scaleHN <- det_curve_sum %>% filter(Distance == 1.5) %>% pull(median) %>% mean() # Approx scaling

y_comb_list <- list()
for(i in 1:max(gs)) {
y_comb_list[[i]] <- colSums(multispecies_data$y[,,i]) %>% 
  as.data.frame() %>%
  `colnames<-`("Count") %>%
  mutate(Distance = midpts,
         Bin = 1:5,
         gs = i,
         Prop = Count/(max(Count)),
         CountS = Count/pa[Bin,i])

}

y_combined <- bind_rows(y_comb_list) %>% 
  group_by(Distance) %>%
  mutate(SumCount = sum(Count)) %>%
  ungroup() %>%
  mutate(propC = Count/SumCount, 
         PropSM = (CountS/max(CountS)) * propC)  %>%
  group_by(Distance) %>% 
  summarise(PropS = sum(PropSM) * scaleHN)

detection_plot_HN <- ggplot(aes(x = Distance), data = det_curve_sum) +
  geom_col(aes(y = PropS), fill = "grey70", colour = "grey20", width = 2.5, data = y_combined, alpha = 0.7) +
  # geom_error bar(aes(ymin = PropSq5, ymax = PropSq95),  data = y_combined) +
  # geom_line(aes(y = HNS)) +
  geom_ribbon(aes(ymin = q5, ymax = q95, fill = herbaceouslvl), alpha = 0.5) +
  geom_line(aes(y = median, colour = herbaceouslvl)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0,1)) +
  scale_fill_manual(values = unname(delwp_cols[c(1,2,3,8)])) +
  scale_colour_manual(values = unname(delwp_cols[c(1,2,3,8)])) +
  ylab("Detection probability (p)") +
  labs(fill = "Herbaceous\nUnderstorey\nCover", colour = "Herbaceous\nUnderstorey\nCover") +
  theme_classic()

detection_plot_HN

Transect-surveys

For Sambar, Fallow and Red Deer we were able to estimate detection rates of the various methods used to detect deer (camera, footprints, pellets, rubbings and wallows).

Show code

all_draws <- model_fits[[top]]$draws("beta_trans_det", format = "df")

# det_marginal_effects <- list()
# det_plot <- list()

obs_vars <- unlist(stringr::str_remove_all(string = labels(terms(transect_formula)), pattern =  "log|scale|\\)|\\("))
obs_cols <- unlist(stringr::str_remove_all(string = colnames(multispecies_data$trans_pred_matrix), pattern =  "log|scale|\\)|\\("))
obs_logs <- unlist(stringr::str_detect(string = colnames(multispecies_data$trans_pred_matrix), pattern =  "log\\("))

obs_labs <- c("Survey Method")

fac <- c(TRUE, TRUE, TRUE, TRUE, TRUE)
fac2 <- c(TRUE)

det_obs <- multispecies_data$transects %>% 
  mutate(Survey = factor(Survey))

det_marginal_effects <- list()
det_plot <- list()

params_w_levs <- levels(det_obs$Survey)

for(i in 1:(length(obs_cols))) {
det_marginal_effects[[i]] <- marginal_effects_cmd(all_draws, 
                                              param = "beta_trans_det", 
                                              param_number = i, log = obs_logs[i],
                                     model_data = det_obs, 
                                     model_column = obs_vars[c(1,attr(multispecies_data$trans_pred_matrix, "assign")[-1])[i]], 
                                     transition = FALSE) %>%
  mutate(group = param, 
         param = params_w_levs[i], 
         factor = fac[i], 
         variable = as.numeric(variable))

if(fac[i]) {
  det_marginal_effects[[i]] <- det_marginal_effects[[i]]
}

}

marginal_prob <- function(x, pwr = 3) {
  xm <- 1-x
  return(1-(xm^pwr))
}

det_marginal_effects_bind <- bind_rows(det_marginal_effects) %>%
  rowwise() %>%
  mutate(value = case_when(param != "Camera" ~ marginal_prob(value), 
                           TRUE ~ value))

det_marginal_effects_split <- split(det_marginal_effects_bind, det_marginal_effects_bind$group)

for(i in 1:length(det_marginal_effects_split)) {
  
  if(fac2[i]) {
    plot_data <- det_marginal_effects_split[[i]] %>% 
      mutate(variable = param, 
             param = group)
  } else {
    plot_data <- det_marginal_effects_split[[i]]
  }

det_plot[[i]] <- plot_data

}
data_for_plot <- bind_rows(det_plot) %>% 
  mutate(species = "All (combined)") %>%
  mutate(Density = "1 deer/km^2^")

data_d3 <- data_for_plot %>% 
  mutate(value = 1 - (1 - value)^3, 
         Density = "3 deer/km^2^")

data_d5 <- data_for_plot %>% 
  mutate(value = 1 - (1 - value)^5, 
         Density = "5 deer/km^2^")

data_combined <- bind_rows(data_for_plot, data_d3, data_d5)

marginal_effects_plot_cmd_all(data_combined, 
                              factor = TRUE,
                              ylab = "Detection probability") +
  xlab("Survey method") +
  scale_y_continuous(limits = c(0,1)) +
  theme(legend.position = "top") +
  scale_fill_manual(values = unname(delwp_cols[c(1,2,3)])) +
  # scale_x_discrete(expand = c(0, 0)) +
  theme_classic() +
  theme(legend.text = ggtext::element_markdown())

Figure 8: Detection Probability for the various methods of survey. Camera trap detection probability is based on average deployment length (53 days), while transects are based on 3 independent transect searches

Abundance Processes

Bioregion

We used a spatially-derived random-effect of the bioregion of the sampled site. This random-effect allows us the ability to make predictions that include the variance associated with this random effect. This random-effect also minimises predictions in areas without deer detection’s (e.g. the mallee).

Show code

bayesplot::color_scheme_set(bayes_theme)
bioregion_contribution <- function(model, data, species, species_index) {
  eps_bioregion <- model$draws("eps_bioregion", format = "matrix")
  
  which_sp <- which(stringr::str_detect(string = colnames(eps_bioregion),
                                         pattern = paste0("eps_bioregion\\[", i)))
  
  eps_bioregion <- eps_bioregion[,which_sp]
  
  bioregion_data <- data[["bioreg_sf"]]
  colnames(eps_bioregion) <- bioregion_data$bioregion
  
  plot <- mcmc_areas(eps_bioregion, prob_outer = 0.9) +
    ggtitle(species)
  
  return(plot)
}

bio_plot_data <- list()

for(i in 1:length(species_names)) {
 bio_plot_data[[i]] <- bioregion_contribution(model = model_fits[[top]], 
                         data = multispecies_data, 
                         species = species_names[i], species_index = i)
} 

cowplot::plot_grid(plotlist = bio_plot_data, ncol = 1)

Figure 9: Influence of the bioregion random effect on abundance (log-scale).

Fixed-parameter Effects

Show code

ab_joined_list <- list()
for(j in 1:length(species_names)) {
beta_draws <- model_fits[[top]]$draws("beta_psi", format = "df")

# which_sp <- which(stringr::str_detect(string = colnames(beta_draws),
#                                     pattern = paste0("beta_psi\\[", i)))
# 
# beta_draws <- beta_draws[,which_sp]

ab_marginal_effects <- list()
ab_plot <- list()

ab_to_use <- ab_formula_3

phi_vars <- unlist(stringr::str_remove_all(string = labels(terms(ab_to_use)), pattern =  "log|scale|sqrt|\\)|\\("))
phi_logs <- unlist(stringr::str_detect(string = labels(terms(ab_to_use)), pattern =  "log\\("))
phi_sqrts <- c(0.5,0.5,1,1,0.5)

phi_labs <- c("Bare Soil (%)", 
              "Nitrogen (%)", 
              "Distance to Pastural Land (m)", 
              "Precipitation Seasonality", 
              'Length of forest edge in 1km2 (m)')

fac <- c(FALSE, FALSE, FALSE, FALSE, FALSE)

for(i in 1:length(phi_vars)) {
ab_marginal_effects[[i]] <- marginal_effects_cmd(beta_draws, 
                                              param = "beta_psi", species_index = j,
                                              param_number = i+1, log = phi_logs[i],
                                     model_data = multispecies_data[["raw_data"]], 
                                     abundance = TRUE,
                                     pwr = phi_sqrts[i],
                                     model_column = phi_vars[i], 
                                     transition = FALSE) %>%
  mutate(species = species_names[j])

}
ab_joined_list[[j]] <- bind_rows(ab_marginal_effects)
}

all_me_data <- bind_rows(ab_joined_list) %>%
  mutate(param = case_when(param == "BareSoil" ~ "Bare Soil (%)", 
                           param == "Nitrogen" ~ "Nitrogen (%)", 
                           param == "PastureDistance" ~ "Distance to Pastural Land (m)", 
                           param == "BIO15" ~ "Precipitation Seasonality", 
                           param == "ForestEdge" ~ 'Amount of Forest Edge per km2 (m)'))
  

marginal_effects_plot_cmd_all(all_me_data, 
                          col = "DarkGreen", 
                          factor = FALSE,
                          ylab = "Contribution to Abundance (log-scale)") +
      ggplot2::facet_grid(species~param, scales = "free") +
  delwp_theme() +
  xlab("Covariate Value")

Figure 10: Marginal response curves for the various fixed-effect parameters used in the model.

Additional Parameters of Interest

Below is a table with several model parameters we estimated:

Show code

key_params <- model_fits[[top]]$summary(c("av_gs", 
                                          "od_mu", 
                                          "od_sd", 
                                          "odRN", 
                                          "activ", 
                                          "Nhat_sum"))

key_params %>% 
  kbl(format = "html", digits = 2, caption = "Several key outputs and parameters of the model not reported earlier") %>%
  kable_styling("striped")

Table 6: Several key outputs and parameters of the model not reported earlier

variable	mean	median	sd	mad	q5	q95	rhat	ess_bulk	ess_tail
av_gs[1]	1.01	1.01	0.00	0.00	1.00	1.02	1.00	1349.70	1322.85
av_gs[2]	1.05	1.04	0.04	0.03	1.01	1.12	1.01	1042.10	1242.36
av_gs[3]	1.09	1.07	0.07	0.05	1.02	1.22	1.00	1976.90	1620.42
av_gs[4]	1.05	1.02	0.07	0.02	1.00	1.17	1.01	592.19	495.72
od_mu[1]	-2.23	-2.23	0.30	0.27	-2.72	-1.75	1.00	1057.37	1220.53
od_mu[2]	-3.49	-3.51	0.39	0.36	-4.11	-2.83	1.01	1063.64	1147.11
od_mu[3]	-1.96	-1.99	0.78	0.76	-3.15	-0.58	1.00	1238.52	1197.40
od_mu[4]	-1.83	-1.96	0.83	0.84	-2.95	-0.29	1.01	776.77	1077.00
od_sd[1]	0.81	0.76	0.30	0.27	0.41	1.37	1.00	813.98	1298.69
od_sd[2]	1.01	0.98	0.44	0.44	0.35	1.82	1.02	352.43	295.37
od_sd[3]	1.67	1.62	0.54	0.54	0.87	2.62	1.00	1476.97	1285.41
od_sd[4]	0.97	0.88	0.62	0.65	0.13	2.08	1.00	890.86	814.25
odRN	0.81	0.79	0.16	0.15	0.59	1.12	1.01	3201.97	1335.53
activ	0.55	0.55	0.02	0.02	0.51	0.59	1.01	3196.18	1106.63
Nhat_sum	270395.70	266181.50	39646.51	37251.81	215378.55	338407.75	1.00	1868.62	1473.08

Relating abundance to vegetation

From our point density estimates of deer, we can model relationships between deer densities at sites and a range of structural vegetation measures collected during the surveys. For this part of the analysis we implement a multivariate regression.

Firstly the response variables we model are:

• Bare Ground Cover (Beta Distribution)
• Native Woody Understorey Cover (Beta Distribution)
• Native Non-woody Understorey Cover (Beta Distribution)
• Presence of seedlings (Negative Binomial Distribution)
• Presence of saplings (Negative Binomial Distribution)
• Presence of exotic flora (Bernoulli Distribution)

While our main predictor of interest is the density of deer at a site, we also control for predictors that are expected to have an impact on vegetation growth and composition. The predictors controlled for include features relating to climate, ecology, topography, fire history and forest structure. They are the following predictors:

• BIO01: Annual Mean Temperature
• BIO04: Temperature Seasonality
  
• BIO06: Min Temperature of Coldest Month
• BIO12: Annual Precipitation
  
• BIO15: Precipitation Seasonality (Coefficient of Variation)
  
• Soil Nitrogen
  
• Time since last fire (factor grouping)
  
• TWIND (Topographic wetness index)
  
• Forest Edge
  
• Canopy Cover
  
• Top Height
• EVC Group (random effect with 19 groups)

Data Preparation

We can use the sum of estimated deer density for the four species at each site as our key predictor variable. Given that we already extracted a range of bioclimatic and environmental covariates for the original model, we can combine those with the data available here.

Note that this model is run on a subset of sites (sites surveyed as part of the statewide deer survey), as the hog deer sites did not have complete vegetation metrics recorded.

Show code

density <- model_fits[[top]]$draws("N_site", format = "matrix")
site_estimates <- list()
site_density <- list()

for(i in 1:n_site) {
  site_estimates[[i]] <- character()
  for(j in 1:length(deer_species_all)) {
    site_estimates[[i]][j] <- paste0("N_site[", j, ",", i, "]")
  }
  site_density[[i]] <- sum(colMeans(density)[site_estimates[[i]]])
}

site_density_ul <- unlist(site_density)

site_density_df <- cams_curated %>%
  mutate(AllDeerDensity = site_density_ul, 
         Bioregion = factor(multispecies_data$site_bioreg), 
         EVC = factor(multispecies_data$site_evc)) %>% 
  filter(ProjectShortName == "StatewideDeer") %>% 
  left_join(multispecies_data[["raw_data"]] %>% select(-EVC, - BIOREGION)) %>%
  mutate(TSLF_bin = as.factor(round(TSLF_bin)))

 site_vars <- dplyr::tbl(con, dbplyr::in_schema("deervic", "curated_site_data")) %>%
    dplyr::filter(SiteID %in% !!c(cams_curated$SiteID, "47191")) %>%
    dplyr::collect() %>%
    dplyr::mutate(HerbaceousUnderstoryCover = NNWHUCover + ENWHUCover,
                  Exotics = EWUCover + ENWHUCover,
                  ExoticPresence = as.integer(as.logical(Exotics)),
                  SeedlingsOrSaplings = as.integer(as.logical(Seedlings+Saplings)),
                  SeedlingsPresence = as.integer(as.logical(Seedlings)),
                  SaplingsPresence = as.integer(as.logical(Saplings)),
                  SiteID = dplyr::case_when(SiteID == "47191" & CameraID == "HO04101053" ~ "47191A",
                                            SiteID == "47191" & CameraID != "HO04101053" ~ "47191B",
                                            TRUE ~ SiteID)) %>% # native + exotic herbaceous cover
    dplyr::arrange(SiteID) %>%
    as.data.frame()
 
 site_density_df_joined <- site_density_df %>% 
   left_join(site_vars) %>%
   mutate(BGroundCover = ifelse(BGroundCover == 0, 0.125, BGroundCover)/100, 
             NWUCover = ifelse(NWUCover == 0, 0.125, NWUCover)/100, 
             NNWHUCover = ifelse(NNWHUCover == 0, 0.125, NNWHUCover)/100)

Model Execution

We use the brms package to run our model in a Bayesian framework (Bürkner 2017, 2018, 2021). Our model is run in parallel on four cores over four chains for 1,000 sampling iterations per chain (1,000 warm up iterations).

Show code

bform1 <- bf(mvbind(BGroundCover, 
                    NWUCover, 
                    NNWHUCover, 
                    Seedlings,
                    Saplings,
                    ExoticPresence) ~ scale(sqrt(AllDeerDensity))
                                    + scale(BIO01)
                                    + scale(BIO04)
                                    + scale(BIO06)
                                    + scale(BIO12)
                                    + scale(BIO15)
                                    + scale(Nitrogen)
                                    + TSLF_bin
                                    + scale(TWIND)
                                    + scale(sqrt(ForestEdge))
                                    + scale(sqrt(CanopyCov))
                                    + scale(sqrt(TopHeight))
                                    + (1|EVC)) 


vegfit <- brm(bform1, 
            data = site_density_df_joined, 
            family = list(Beta(link = "logit"), 
                          Beta(link = "logit"),
                          Beta(link = "logit"),
                          negbinomial(),
                          negbinomial(),
                          bernoulli(link = "logit")),
            inits = "0",
            chains = 4, 
            cores = 4, 
            iter = 2000, 
            warmup = 1000,
            backend = "cmdstanr")

Model Results

Model Summary

Below we show a complete summary of the model, including all coefficients estimated for all the parameters and response variables. Standard deviations for the EVC random effect and the estimates of the \(\phi\) parameter in the beta-distribution models, and the shape parameter in the negative binomial models are also shown. Note that all Rhats > 1.05 and all ESS > 200, suggesting good chain mixing and convergence.

Show code

summary(vegfit) 

 Family: MV(beta, beta, beta, negbinomial, negbinomial, bernoulli) 
  Links: mu = logit; phi = identity
         mu = logit; phi = identity
         mu = logit; phi = identity
         mu = log; shape = identity
         mu = log; shape = identity
         mu = logit 
Formula: BGroundCover ~ scale(sqrt(AllDeerDensity)) + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
         NWUCover ~ scale(sqrt(AllDeerDensity)) + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
         NNWHUCover ~ scale(sqrt(AllDeerDensity)) + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
         Seedlings ~ scale(sqrt(AllDeerDensity)) + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
         Saplings ~ scale(sqrt(AllDeerDensity)) + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
         ExoticPresence ~ scale(sqrt(AllDeerDensity)) + scale(BIO01) + scale(BIO04) + scale(BIO06) + scale(BIO12) + scale(BIO15) + scale(Nitrogen) + TSLF_bin + scale(TWIND) + scale(sqrt(ForestEdge)) + scale(sqrt(CanopyCov)) + scale(sqrt(TopHeight)) + (1 | EVC) 
   Data: site_density_df_joined (Number of observations: 252) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~EVC (Number of levels: 19) 
                             Estimate Est.Error l-95% CI u-95% CI
sd(BGroundCover_Intercept)       0.27      0.11     0.06     0.52
sd(NWUCover_Intercept)           0.13      0.09     0.01     0.34
sd(NNWHUCover_Intercept)         0.24      0.13     0.02     0.53
sd(Seedlings_Intercept)          1.24      0.42     0.54     2.22
sd(Saplings_Intercept)           0.86      0.32     0.25     1.55
sd(ExoticPresence_Intercept)     0.54      0.37     0.03     1.40
                             Rhat Bulk_ESS Tail_ESS
sd(BGroundCover_Intercept)   1.00     1489     1529
sd(NWUCover_Intercept)       1.00     1658     1929
sd(NNWHUCover_Intercept)     1.00      988     1110
sd(Seedlings_Intercept)      1.00      952     1388
sd(Saplings_Intercept)       1.00      996     1036
sd(ExoticPresence_Intercept) 1.01      956     1851

Population-Level Effects: 
                                       Estimate Est.Error l-95% CI
BGroundCover_Intercept                    -2.23      0.21    -2.65
NWUCover_Intercept                        -0.64      0.17    -0.99
NNWHUCover_Intercept                      -0.30      0.19    -0.69
Seedlings_Intercept                        3.34      0.53     2.30
Saplings_Intercept                         3.47      0.45     2.59
ExoticPresence_Intercept                  -1.37      0.50    -2.35
BGroundCover_scalesqrtAllDeerDensity      -0.19      0.09    -0.37
BGroundCover_scaleBIO01                    0.35      0.42    -0.46
BGroundCover_scaleBIO04                    0.08      0.26    -0.42
BGroundCover_scaleBIO06                   -0.22      0.48    -1.14
BGroundCover_scaleBIO12                   -0.20      0.15    -0.50
BGroundCover_scaleBIO15                   -0.22      0.08    -0.38
BGroundCover_scaleNitrogen                -0.03      0.16    -0.35
BGroundCover_TSLF_bin3                     0.13      0.22    -0.28
BGroundCover_TSLF_bin4                     0.09      0.24    -0.38
BGroundCover_TSLF_bin5                    -0.21      0.25    -0.69
BGroundCover_scaleTWIND                    0.05      0.11    -0.16
BGroundCover_scalesqrtForestEdge          -0.01      0.07    -0.15
BGroundCover_scalesqrtCanopyCov           -0.01      0.07    -0.15
BGroundCover_scalesqrtTopHeight           -0.25      0.08    -0.40
NWUCover_scalesqrtAllDeerDensity          -0.17      0.08    -0.33
NWUCover_scaleBIO01                        0.03      0.42    -0.82
NWUCover_scaleBIO04                       -0.02      0.26    -0.54
NWUCover_scaleBIO06                       -0.03      0.48    -0.98
NWUCover_scaleBIO12                        0.22      0.15    -0.07
NWUCover_scaleBIO15                        0.18      0.08     0.02
NWUCover_scaleNitrogen                    -0.02      0.16    -0.33
NWUCover_TSLF_bin3                        -0.38      0.19    -0.76
NWUCover_TSLF_bin4                        -0.59      0.23    -1.03
NWUCover_TSLF_bin5                        -0.93      0.25    -1.40
NWUCover_scaleTWIND                       -0.04      0.12    -0.26
NWUCover_scalesqrtForestEdge               0.06      0.08    -0.09
NWUCover_scalesqrtCanopyCov               -0.05      0.07    -0.19
NWUCover_scalesqrtTopHeight               -0.13      0.08    -0.28
NNWHUCover_scalesqrtAllDeerDensity         0.22      0.09     0.06
NNWHUCover_scaleBIO01                      0.28      0.49    -0.73
NNWHUCover_scaleBIO04                     -0.35      0.30    -0.92
NNWHUCover_scaleBIO06                     -0.58      0.56    -1.64
NNWHUCover_scaleBIO12                      0.37      0.16     0.05
NNWHUCover_scaleBIO15                      0.12      0.09    -0.05
NNWHUCover_scaleNitrogen                  -0.27      0.17    -0.60
NNWHUCover_TSLF_bin3                      -0.29      0.21    -0.70
NNWHUCover_TSLF_bin4                      -0.14      0.24    -0.61
NNWHUCover_TSLF_bin5                      -0.15      0.25    -0.63
NNWHUCover_scaleTWIND                     -0.08      0.13    -0.33
NNWHUCover_scalesqrtForestEdge             0.00      0.08    -0.16
NNWHUCover_scalesqrtCanopyCov              0.03      0.08    -0.12
NNWHUCover_scalesqrtTopHeight              0.02      0.08    -0.15
Seedlings_scalesqrtAllDeerDensity          0.27      0.22    -0.16
Seedlings_scaleBIO01                       0.07      0.93    -1.73
Seedlings_scaleBIO04                      -0.08      0.59    -1.27
Seedlings_scaleBIO06                      -0.47      1.02    -2.48
Seedlings_scaleBIO12                       0.56      0.45    -0.31
Seedlings_scaleBIO15                      -0.16      0.29    -0.73
Seedlings_scaleNitrogen                   -0.93      0.40    -1.71
Seedlings_TSLF_bin3                       -1.15      0.52    -2.13
Seedlings_TSLF_bin4                       -1.19      0.53    -2.24
Seedlings_TSLF_bin5                       -1.61      0.58    -2.72
Seedlings_scaleTWIND                      -0.10      0.25    -0.59
Seedlings_scalesqrtForestEdge             -0.39      0.17    -0.73
Seedlings_scalesqrtCanopyCov              -0.09      0.18    -0.43
Seedlings_scalesqrtTopHeight              -0.93      0.20    -1.34
Saplings_scalesqrtAllDeerDensity          -0.07      0.17    -0.41
Saplings_scaleBIO01                        0.26      0.81    -1.28
Saplings_scaleBIO04                        0.10      0.50    -0.87
Saplings_scaleBIO06                        0.17      0.89    -1.54
Saplings_scaleBIO12                        0.31      0.32    -0.31
Saplings_scaleBIO15                       -0.44      0.19    -0.81
Saplings_scaleNitrogen                     0.13      0.33    -0.51
Saplings_TSLF_bin3                        -0.57      0.46    -1.46
Saplings_TSLF_bin4                        -1.75      0.48    -2.68
Saplings_TSLF_bin5                        -1.31      0.54    -2.35
Saplings_scaleTWIND                        0.09      0.23    -0.35
Saplings_scalesqrtForestEdge              -0.30      0.15    -0.60
Saplings_scalesqrtCanopyCov                0.02      0.14    -0.27
Saplings_scalesqrtTopHeight               -0.56      0.17    -0.90
ExoticPresence_scalesqrtAllDeerDensity     0.60      0.22     0.18
ExoticPresence_scaleBIO01                  0.36      1.11    -1.77
ExoticPresence_scaleBIO04                 -0.06      0.71    -1.47
ExoticPresence_scaleBIO06                 -0.60      1.31    -3.24
ExoticPresence_scaleBIO12                  0.23      0.40    -0.52
ExoticPresence_scaleBIO15                  0.38      0.21    -0.03
ExoticPresence_scaleNitrogen              -1.21      0.46    -2.12
ExoticPresence_TSLF_bin3                   0.68      0.57    -0.44
ExoticPresence_TSLF_bin4                   0.70      0.64    -0.59
ExoticPresence_TSLF_bin5                   1.00      0.64    -0.32
ExoticPresence_scaleTWIND                 -0.40      0.30    -0.99
ExoticPresence_scalesqrtForestEdge         0.08      0.19    -0.29
ExoticPresence_scalesqrtCanopyCov          0.24      0.19    -0.13
ExoticPresence_scalesqrtTopHeight          0.09      0.21    -0.31
                                       u-95% CI Rhat Bulk_ESS
BGroundCover_Intercept                    -1.84 1.00     2168
NWUCover_Intercept                        -0.31 1.00     2814
NNWHUCover_Intercept                       0.07 1.00     2853
Seedlings_Intercept                        4.38 1.00     1633
Saplings_Intercept                         4.33 1.00     1895
ExoticPresence_Intercept                  -0.38 1.00     2673
BGroundCover_scalesqrtAllDeerDensity      -0.03 1.00     4322
BGroundCover_scaleBIO01                    1.16 1.00     1904
BGroundCover_scaleBIO04                    0.60 1.00     1988
BGroundCover_scaleBIO06                    0.69 1.00     1939
BGroundCover_scaleBIO12                    0.09 1.00     4321
BGroundCover_scaleBIO15                   -0.06 1.00     3184
BGroundCover_scaleNitrogen                 0.29 1.00     3855
BGroundCover_TSLF_bin3                     0.56 1.00     2314
BGroundCover_TSLF_bin4                     0.55 1.00     2347
BGroundCover_TSLF_bin5                     0.29 1.00     2378
BGroundCover_scaleTWIND                    0.27 1.00     5203
BGroundCover_scalesqrtForestEdge           0.14 1.00     5038
BGroundCover_scalesqrtCanopyCov            0.14 1.00     5102
BGroundCover_scalesqrtTopHeight           -0.10 1.00     5962
NWUCover_scalesqrtAllDeerDensity          -0.03 1.00     5148
NWUCover_scaleBIO01                        0.84 1.00     2005
NWUCover_scaleBIO04                        0.50 1.00     2031
NWUCover_scaleBIO06                        0.94 1.00     1966
NWUCover_scaleBIO12                        0.51 1.00     4152
NWUCover_scaleBIO15                        0.34 1.00     3601
NWUCover_scaleNitrogen                     0.29 1.00     4203
NWUCover_TSLF_bin3                         0.01 1.00     3003
NWUCover_TSLF_bin4                        -0.13 1.00     2853
NWUCover_TSLF_bin5                        -0.43 1.00     2971
NWUCover_scaleTWIND                        0.20 1.00     4799
NWUCover_scalesqrtForestEdge               0.21 1.00     5152
NWUCover_scalesqrtCanopyCov                0.10 1.00     5259
NWUCover_scalesqrtTopHeight                0.02 1.00     5646
NNWHUCover_scalesqrtAllDeerDensity         0.39 1.00     4935
NNWHUCover_scaleBIO01                      1.19 1.00     1551
NNWHUCover_scaleBIO04                      0.27 1.00     1584
NNWHUCover_scaleBIO06                      0.53 1.00     1560
NNWHUCover_scaleBIO12                      0.70 1.00     3475
NNWHUCover_scaleBIO15                      0.30 1.00     3666
NNWHUCover_scaleNitrogen                   0.07 1.00     3849
NNWHUCover_TSLF_bin3                       0.15 1.00     2711
NNWHUCover_TSLF_bin4                       0.33 1.00     2935
NNWHUCover_TSLF_bin5                       0.36 1.00     2596
NNWHUCover_scaleTWIND                      0.17 1.00     3906
NNWHUCover_scalesqrtForestEdge             0.16 1.00     5144
NNWHUCover_scalesqrtCanopyCov              0.18 1.00     4698
NNWHUCover_scalesqrtTopHeight              0.18 1.00     3936
Seedlings_scalesqrtAllDeerDensity          0.71 1.00     3494
Seedlings_scaleBIO01                       1.90 1.00     1565
Seedlings_scaleBIO04                       1.04 1.00     1746
Seedlings_scaleBIO06                       1.48 1.00     1595
Seedlings_scaleBIO12                       1.49 1.00     1904
Seedlings_scaleBIO15                       0.43 1.00     2270
Seedlings_scaleNitrogen                   -0.14 1.00     4452
Seedlings_TSLF_bin3                       -0.13 1.00     2474
Seedlings_TSLF_bin4                       -0.15 1.00     2221
Seedlings_TSLF_bin5                       -0.47 1.00     2224
Seedlings_scaleTWIND                       0.39 1.00     3242
Seedlings_scalesqrtForestEdge             -0.05 1.00     4653
Seedlings_scalesqrtCanopyCov               0.28 1.00     4454
Seedlings_scalesqrtTopHeight              -0.53 1.00     4612
Saplings_scalesqrtAllDeerDensity           0.26 1.00     4647
Saplings_scaleBIO01                        1.84 1.00     1365
Saplings_scaleBIO04                        1.07 1.00     1496
Saplings_scaleBIO06                        1.86 1.00     1406
Saplings_scaleBIO12                        0.95 1.00     2389
Saplings_scaleBIO15                       -0.06 1.00     3288
Saplings_scaleNitrogen                     0.80 1.00     3575
Saplings_TSLF_bin3                         0.34 1.00     2091
Saplings_TSLF_bin4                        -0.80 1.00     2307
Saplings_TSLF_bin5                        -0.27 1.00     2216
Saplings_scaleTWIND                        0.55 1.00     3646
Saplings_scalesqrtForestEdge              -0.01 1.00     5397
Saplings_scalesqrtCanopyCov                0.30 1.00     4268
Saplings_scalesqrtTopHeight               -0.23 1.00     3460
ExoticPresence_scalesqrtAllDeerDensity     1.04 1.00     4137
ExoticPresence_scaleBIO01                  2.66 1.00     2060
ExoticPresence_scaleBIO04                  1.33 1.00     2014
ExoticPresence_scaleBIO06                  1.96 1.00     1978
ExoticPresence_scaleBIO12                  1.04 1.00     4530
ExoticPresence_scaleBIO15                  0.80 1.00     3909
ExoticPresence_scaleNitrogen              -0.34 1.00     4093
ExoticPresence_TSLF_bin3                   1.78 1.00     2519
ExoticPresence_TSLF_bin4                   1.90 1.00     2431
ExoticPresence_TSLF_bin5                   2.23 1.00     2680
ExoticPresence_scaleTWIND                  0.20 1.00     5473
ExoticPresence_scalesqrtForestEdge         0.46 1.00     5535
ExoticPresence_scalesqrtCanopyCov          0.62 1.00     5583
ExoticPresence_scalesqrtTopHeight          0.48 1.00     5344
                                       Tail_ESS
BGroundCover_Intercept                     2571
NWUCover_Intercept                         2971
NNWHUCover_Intercept                       2725
Seedlings_Intercept                        2330
Saplings_Intercept                         2639
ExoticPresence_Intercept                   2984
BGroundCover_scalesqrtAllDeerDensity       3086
BGroundCover_scaleBIO01                    2511
BGroundCover_scaleBIO04                    2622
BGroundCover_scaleBIO06                    2585
BGroundCover_scaleBIO12                    3360
BGroundCover_scaleBIO15                    3109
BGroundCover_scaleNitrogen                 3158
BGroundCover_TSLF_bin3                     2833
BGroundCover_TSLF_bin4                     2942
BGroundCover_TSLF_bin5                     2915
BGroundCover_scaleTWIND                    2841
BGroundCover_scalesqrtForestEdge           3392
BGroundCover_scalesqrtCanopyCov            3387
BGroundCover_scalesqrtTopHeight            3337
NWUCover_scalesqrtAllDeerDensity           3118
NWUCover_scaleBIO01                        2275
NWUCover_scaleBIO04                        2380
NWUCover_scaleBIO06                        2356
NWUCover_scaleBIO12                        3185
NWUCover_scaleBIO15                        3059
NWUCover_scaleNitrogen                     3217
NWUCover_TSLF_bin3                         2968
NWUCover_TSLF_bin4                         2969
NWUCover_TSLF_bin5                         2932
NWUCover_scaleTWIND                        3447
NWUCover_scalesqrtForestEdge               3293
NWUCover_scalesqrtCanopyCov                3107
NWUCover_scalesqrtTopHeight                2509
NNWHUCover_scalesqrtAllDeerDensity         2983
NNWHUCover_scaleBIO01                      2217
NNWHUCover_scaleBIO04                      2313
NNWHUCover_scaleBIO06                      2226
NNWHUCover_scaleBIO12                      2976
NNWHUCover_scaleBIO15                      3302
NNWHUCover_scaleNitrogen                   3502
NNWHUCover_TSLF_bin3                       2810
NNWHUCover_TSLF_bin4                       2975
NNWHUCover_TSLF_bin5                       2449
NNWHUCover_scaleTWIND                      3123
NNWHUCover_scalesqrtForestEdge             3358
NNWHUCover_scalesqrtCanopyCov              2970
NNWHUCover_scalesqrtTopHeight              2458
Seedlings_scalesqrtAllDeerDensity          2744
Seedlings_scaleBIO01                       2334
Seedlings_scaleBIO04                       2462
Seedlings_scaleBIO06                       2237
Seedlings_scaleBIO12                       2162
Seedlings_scaleBIO15                       2552
Seedlings_scaleNitrogen                    2978
Seedlings_TSLF_bin3                        2746
Seedlings_TSLF_bin4                        2817
Seedlings_TSLF_bin5                        2774
Seedlings_scaleTWIND                       2904
Seedlings_scalesqrtForestEdge              2959
Seedlings_scalesqrtCanopyCov               3287
Seedlings_scalesqrtTopHeight               2983
Saplings_scalesqrtAllDeerDensity           3516
Saplings_scaleBIO01                        2091
Saplings_scaleBIO04                        2088
Saplings_scaleBIO06                        2289
Saplings_scaleBIO12                        2628
Saplings_scaleBIO15                        3319
Saplings_scaleNitrogen                     2775
Saplings_TSLF_bin3                         2918
Saplings_TSLF_bin4                         2919
Saplings_TSLF_bin5                         3048
Saplings_scaleTWIND                        3253
Saplings_scalesqrtForestEdge               3316
Saplings_scalesqrtCanopyCov                3161
Saplings_scalesqrtTopHeight                2901
ExoticPresence_scalesqrtAllDeerDensity     3079
ExoticPresence_scaleBIO01                  2150
ExoticPresence_scaleBIO04                  1661
ExoticPresence_scaleBIO06                  1980
ExoticPresence_scaleBIO12                  3115
ExoticPresence_scaleBIO15                  3066
ExoticPresence_scaleNitrogen               3246
ExoticPresence_TSLF_bin3                   3246
ExoticPresence_TSLF_bin4                   3013
ExoticPresence_TSLF_bin5                   3140
ExoticPresence_scaleTWIND                  2940
ExoticPresence_scalesqrtForestEdge         2887
ExoticPresence_scalesqrtCanopyCov          3097
ExoticPresence_scalesqrtTopHeight          2736

Family Specific Parameters: 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
phi_BGroundCover     6.89      0.73     5.57     8.45 1.00     4735
phi_NWUCover         4.10      0.37     3.42     4.87 1.00     6595
phi_NNWHUCover       3.08      0.26     2.60     3.61 1.00     6327
shape_Seedlings      0.29      0.04     0.23     0.37 1.00     2753
shape_Saplings       0.38      0.04     0.30     0.47 1.00     3974
                 Tail_ESS
phi_BGroundCover     3057
phi_NWUCover         2461
phi_NNWHUCover       3052
shape_Seedlings      2862
shape_Saplings       2721

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Deer Density Conditional Effects

Given our main interest in the model is deer density and how it impacts the suite of response variables, we generate marginal/conditional effects plots for the relationship between deer and the vegetation metric. Based on these plots we can see a “significant” (confidence intervals not-overlapping zero) effect of deer density on: bare ground (-ve), native woody understorey (-ve), native non-woody/herbaceous understorey (+ve), and presence of exotic species (+ve).

Show code

plot_list <- conditional_effects(vegfit, effects = "AllDeerDensity", prob = 0.9)
plot_list_50 <- conditional_effects(vegfit, effects = "AllDeerDensity", prob = 0.5)

plot_deer_relationship <- function(plot_data, plot_data_50, yaxis) {
  ggplot(plot_data) +
    geom_ribbon(aes(x = AllDeerDensity, ymin = `lower__`, ymax = `upper__`), 
                fill = delwp_cols[["Navy"]], alpha = 0.3) +
    geom_ribbon(aes(x = AllDeerDensity, ymin = `lower__`, ymax = `upper__`), 
                fill = delwp_cols[["Navy"]], alpha = 0.3, data = plot_data_50) +
    geom_smooth(aes(x = AllDeerDensity, y = `estimate__`), 
                colour = delwp_cols[["Navy"]], method = "loess", formula = y~x) +
    ylab(paste(yaxis)) +
    delwp_theme()
    
}

cowplot::plot_grid(plot_deer_relationship(plot_list[[1]], plot_list_50[[1]], "Bare Ground Cover"), 
                   plot_deer_relationship(plot_list[[2]], plot_list_50[[2]], "Native Woody Understorey Cover"),
                   plot_deer_relationship(plot_list[[3]], plot_list_50[[3]], "Native Herbaceous Understorey Cover"),
                   plot_deer_relationship(plot_list[[4]], plot_list_50[[4]], "Seedling Counts") + scale_y_sqrt(),
                   plot_deer_relationship(plot_list[[5]], plot_list_50[[5]], "Sapling Counts"),
                   plot_deer_relationship(plot_list[[6]], plot_list_50[[6]], "Presence of Exotic Flora"), 
                   ncol = 2, labels = "AUTO")

Figure 11: Impact of deer density on a range of vegetation measure, including (A) bare ground cover, (B) native woody understorey, (C) native herbaceous understorey cover, (D) seedlings, (E) saplings, and (F) th presence of exotic species. The line shows the mean estimate, with the shading representing 50% and 90% confidence intervals.